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In this paper bounds for the error of a computed inverse of a matrix are developed. These are 
then applied to the solution of a single system. Methods for improving the approximate inverse are 
then discussed with some observations on the dangers involved in their practical use on a computer 
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Introduction 

The advent of modern high-speed electronic computers has led to the solution of many prob- 
lems that had been considered almost unsolvable some years ago. The essential contribution of 
these computers is that they can do simple things extremely fast. Mathematical problems, like 
the solution of a system of many linear algebraic equations, can be done in a few seconds whereas 
it would take many hours by a person with a desk calculator and even longer if he used only pencil 
and paper. If we actually solved the problem in either of the latter ways, we would have a good 
idea of the errors we allowed to enter the process and could keep track of them. It is in this area 
that high-speed computers have some limitations. We feed numbers into a computer and the com- 
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puter prints out other numbers — supposedly the solution of our problem. But how do we know that 
these numbers are the solution of our problem? That this is a problem is pointed out very well by 
Leslie Fox: 

In the desk-machine era of Comrie and Peters, the great compilers of mathematical tables, the possibility of errors 
and the perpetual need for guarding against them were fully realised. They would be appalled by the self-confidence 
of their successors! [1, p. 16] 1 

It certainly is possible that the machine could malfunction; some physical part could be defective — 
this has, in fact, happened and with no indication to the operator. Some modern text books on 
numerical analysis do not seem to recognize this possibility and this is most unfortunate [2, p. 111]. 
We will indicate in section 6 one way of checking for this kind of error. 

Another possibility that exists is that our instructions to the machine (the "program") may 
not accomplish exactly what we had intended. There are many reasons for this, such as human 
error in punching the deck, the rounding process that continually occurs, subtraction of numbers 
that are very close to each other, and even poor programming. These difficulties are being investi- 
gated very seriously now by many men in the field, especially some in the Special Interest Group 
in Numerical Mathematics of the Association for Computing Machinery. We must be very sure 
that the computer program we are using does, in fact, perform successfully. Evaluation of existing 
computer programs is going on [3] and, in fact, we have devoted the second part of our research 
to this very thing, but, of course, only in a specific area. The criteria for this evaluation will be 
discussed at the beginning of section 6 and the search for these criteria gave the impetus for the 
first part of this research. 

In this first part we derive some theoretical bounds for the accuracy of a computed solution 
of a system of linear algebraic equations. We concentrate first on the problem of finding an error 
bound for the computed inverse of a matrix (the solution of a set of linear algebraic systems) and 
then consider the solution of a single system. (This order may seem backwards at first, but we indi- 
cate its necessity at the end of section 2.) Our orientation will be to use what is computationally 
available. That is to say, if we are interested in the accuracy of a computed inverse, X, we will 
assume that the residual, Y=I — AX, is also available. We will then discuss ways of improving the 
computed inverse and the calculated solution of a single system and finally discuss error bounds 
for some special matrices. 

In spite of the power of a high-speed computer, these theoretical relations discussed in the 
first part are of the utmost importance. We attempt to derive realistic (practical) bounds and we 
have opportunity in section 3 to show where the computed results must be very carefully examined; 
that is, they contradict our theoretical relations. Hence, our theoretical results can always be used 
at least in this negative way. 

One of the concepts basic to a theoretical treatment of errors is that of a norm, and we feel it 
would be useful to give a brief discussion of this concept here, before proceeding to the body of 
the text. A full treatment of this topic can be found in Computational Methods of Linear Algebra, by 
Faddeev and Faddeeva. In dealing with vectors and matrices we frequently need some way of 
describing the concept of the size of a vector, x, or a matrix, A. We shall use the notation N(x) or 
N(A) to indicate this concept. 

In two-dimensional Euclidean space, the size of a vector* = (xi, #2) is just the length (always 
non-negative) of that vector: 

N(x) = {x\ + x\)V 2 . 

The length of such a vector has some immediate properties: 

1. N{x) > for x # and 7V(0) = 0; 

2. N(cx) = \c\N(x) for any numerical multiplier c; 

3. N(x + y) *N(x) +N(y). 



1 Figures in brackets indicate the literature references on p. 309. 
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The usual length of a vector is just one way of describing its size. Other ways would depend upon 
our main point of emphasis. If we are interested in the distance in reference to the coordinate axes, 
we might take the norm of x to be the maximum distance from the origin along one of the axis: 

N(x) = max \xi\, 

i 

For this definition, the same three properties also hold. The main thing, then, that we wish to 
associate with the concept of a norm is that it possess the three properties mentioned above. The 
generalization to n dimensions of each of the two examples given is obvious (although the proofs of 
some of the properties may not be), so let us consider the norm of a matrix. 

Definition: The norm of a square matrix A is a nonnegative number N(A) which satisfies the 
conditions 

1. N(A) >0,i/A^0 and N(0) = 0; 

2. N(cA)= |c|N(A); 

3. N(A + B) *£ N(A) + N(B) 

4. N(AB) ^ N(A)N(B). 

Some examples are in order. Similar to the length of a vector, we could take the square root of the 
sum of the squares of all the elements of A ; 



N{A) = f ^ \ofj\ 2 I (the Frobenius norm). 



Corresponding to the maximum element of a vector, we could take the maximum row sum of A: 

N(A) = max]£|ay|. 

' j 

Another one frequently used is also related to the maximum element of an n X n matrix: 

N(A) = n max |ay|. 

' . j 

In this case, the multiplier n is needed in order to satisfy the required properties. 

Property 4 is peculiar to matrices and we have one special case to consider here; viz., the 
product Ax, a matrix times a vector. In this case we must choose the norms to be used in such a 
way that 

N(Ax) ^N(A)N(x). 

This is called the compatibility condition. For each norm of a vector, there may be more than one 
norm of a matrix that can be used for N(A) but not every matrix norm may satisfy this inequality. 
Our main need for these norms arises when we consider the error of an approximate inverse 
X. We need some way to discuss the "size" of the error, A -1 — X , and the concept of a norm does 
this very well. 

1. Error Bounds for the Approximate Inverse 

In working with norms, frequent use if made of the property that N(AB) ^N(A)N(B) . Before 
discussing various error bounds, it will be useful to examine this inequality, as it is the main reason 
why so many error estimates are overestimates. Let us take the Frobenius norm: 



A^) = (2 k;l 2 ) 



1/2 
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Then, for a general A and B, 

N(AB)* = "2 %a ik b kj 



^21 2 l a *^l 

i,j V ft 
* > j > r , s 

=s Y -(|a <r 6„|» + |a lf M*) 

i ,j , r , s ^ 

= 2 l a *> 6 «"| 2 
i> j , r , s 

=(l\^)h\b sj \A 

i, r s,j ' 

=N(A)*N(B) 2 . 

From this derivation we can determine when equality is possible. If we consider matrices with 
only non-negative elements, the first inequality in the above becomes an equality. Hence, the 
remaining inequality results from the fact that 2ab ^ a 2 + b 2 for all a and b. This is an equality if 
and only if a= b. In our problem, this means if and only if ai r b S j = a\ s b r j or if and only if 

air = b r j 
cits b sj 

which says that the columns of A are in the same ratio as the rows of B. This conclusion can be 
written as b r j — (b S jlai S )ai r for all r; that is, each column of B is a multiple of a row of A. The last 
statement can be seen very clearly if we examine the Cauchy-Bunyakovskii inequality which is the 
basis for the step in the above derivation that we are discussing. This inequality states that for 
two vectors, X and Y, we have | (Z, Y) | ^ \X\- \Y\ where (X, Y) is the scalar product. The familiar 
proof is as follows. Let us consider a nonzero X, and let Z = Y — aX, where a = (Y, X)f (X, X) . We 
note that (Z, X) = (F, X) - a(Z, X) = 0. Then 

|Z|2 = (Z,Z)= (Z,Y-aX)= (Z, Y) = (Y - aX , Y) = (Y,Y)-a(X,Y) 

_ (Y,X)(X,Y) JX\*\Y\*-\(X,Y)\* 
1 ' } (X 9 X) \x\* 

Consequently |Z| 2 |F| 2 — | (X, Y) | 2 = |Z| 2 |Z| 2 ^ 0, and the desired inequality follows. We will have 
equality if and only if Z = 0. This says that Y=aX; i.e., one vector is a scalar multiple of the other. 

Although this knowledge of when equality occurs does not give us any idea of the magnitude 
of the inequality, it certainly indicates that N(AB) will seldom be equal to N(A)N(B). Since we 
have no knowledge of when the elements of A and B are going to be non-negative, the first inequality 
in the derivation just makes matters worse. If the signs of the elements of A and B are such that 
cancellation takes place, this inequality could also cause a tremendous difference. It has been our 
experience in general that the norm of a product is very much less than the product of the norms. 
(See the example at the end of this section.) 

There is one case, however, when equality will always hold if yV is the Frobenius norm. If U 
and V are unitary matrices, then N(UAV) = N(A). To see this, we use N(A) 2 = trace {AA*), A* 
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being the conjugate transpose of A, and the fact that if two matrices are similar, their traces are 
equal (since the trace is equal to the sum of the eigenvalues). Thus 

N(UAV)- = tr(UAVV*A*U*) 
= tr(UAA*U*) 
= tr(AA*) 

= N{A)\ 

This equality has the following two consequences. If A is a normal matrix, it is unitarily equiva- 
lent to a diagonal matrix, and we have N(AB) = N(U*DUB) = N(DUB) ^ N(D)N(UB) 
= N(D)N(B) ^ n ]l ' 2 max \i,\N(B) where It are the eigenvalues of A. The second one also concerns 
the eigenvalues of A. We use a result of Schur which states that an arbitrary complex matrix A 
can be transformed into triangular form by a unitary matrix U. Now for any triangular matrix 7 1 , 
we have 

X M* * N(Ty 

i 

with equality if and only if T is diagonal. Hence, 

2 \h\ 2 = 2 \tu\ 2 « N(Ty = N(U*AUY = N{AY 

j i 

with equality if and only if A is unitarily equivalent to a diagonal matrix; i.e., when A is normal. 

As a final note on the relation between N(AB) and N(A)N(B), we would like to mention that 
there exists no norm such that equality holds for all A and B. To prove this, we would just consider 
an A = B such that A ¥> but A- = 0. 

Our main concern will be with bounds furN(A~ l ) and N(A X ~X) where X is an approximate 
inverse and /V is any norm. To derive upper bounds for N(A l ) is a nontrivial matter. Since A may 
be arbitrarily close to a singular matrix, uniform upper bounds for iVC^"' 1 ) are not to be found. 
Hence, we must make use of information that can be readily available. The most useful quantity 
for an approximate inverse X is the residual matrix Y = I — AX. The "size" of F, of course, depends 
on how good an approximation X is but we have found that N(Y) is invariably less than 1 and, in 
fact, actually less than J. 

The fact that N(Y) is computationally less than 1 turns out to be quite important because this 
inequality is basic to practically all our theoretical results. Its main use follows from the following 
two lemmas. 

Lemma 1.1. If B is a matrix such that N(B) < 1, then 

(a) (I ± B) is nonsingular and 

(b) N[(I±B)->] «^JL_(fN(D = l. 

PROOF: (a) Assume {I±B) is a singular. Then there exists some vector x t^ such that 

(I±B)x = 
or 

x = TBx 
and 

N(x) ^N(B)N{x). 

Since N(x) > 0, this last inequality means that N(B) ^ 1 which contradicts our condition that 
N(B) < 1. 
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(b) Since (7±B) is nonsingular, we can let 

R=(I±B)-K 



Then 



and 



R- l =(I±B) 
I = R±RB 
R=I+RB 
N(R) ^N(I)+N(R)N(B) 
Then *1+N(R)N(B)XN(I) = 1. 

N(R)[l-N(B)]^l 

1 



and 



N(R) 



l-N(B) 



Another relation that needs the same condition and one that is also very useful is the following. 
Lemma 1.2. If B is a matrix such that N(B) < 1, then 

(I-B) , = I + B-hB 2 + . . . +B n + .... 

Proof: Let P n = (/+ B + S 2 + . . . +B"- l )(I-B) 

= I-B». 
Then I-P n = B» 

and N(I-P tl )=N(B») 

^N(B)". 

Since N(B) < 1, N(B) n approaches zero as n goes to infinity so that N(I — P n ) approaches zero 
also and I — P n must converge to the matrix, which means P n converges to the identity matrix 
and hence (I + B + B 2 + . . . +B"~ l ) converges to (I-B)\ 

With these two lemmas as background, we can proceed to derive the bounds for the matrices 
in which we are interested. 

THEOREM 1.1. Let X be an approximate inverse of A and let Y = I — AX. Then i/N(Y) < 1, 

N(X) gNfA-X N(X) • 
1 + N(Y) y ' l-N(Y) 

Proof: We shall orove the upper bound first. 

Y = I - AX 
AX = I-Y 

X~ t A- , = (I-Y)~ l 

A- ] =X(I-Y)-> 

and aky\ 

K(A-<)^ l _ ( *> Y) i{N(I) = l. 

We may remove the restriction on N{I) by the following: 

A- l =X(I-Y)- 1 

= X(I + Y+Y* + Y 3 + . . .) 
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= X + XY + XY 1 + XY : < + . . . 
and N(A-i) ^ N(X) +N(X)N(Y) +N(X)N(T i ) + . . . 

as N(X) +N(X)N(Y) +N(X)N(Y)-+ . . . 

= N(X)[] +N(Y) +N(Y)-+ . . .] 

N(X) 
\-N{Y) 

For the lower bound, we proceed in a similar way. 

Y = I - AX 

AX = I-Y 

X = A~'(I-Y) 

= A-'-A~ f Y 

and N(X) ^ N(A-*)(1+N(Y)) 

so that „t/-«a 

l + N(Y) { } 

We note that taking Y = I — XA, the lefthand residual matrix, would yield identical results. 
COROLLARY. For any nonsingular A, let D be the diagonal and let — M represent the off- 
diagonal part of A. 
Then, i/N(D" 1 M)<l, 

N (D - 1 ) ^ N(A-) ^ N ( D "> ■ 
1 + N(D" 1 M) v ; l-N(D-'M) 

Proof: If we write A=D-M = D(I-D~ l M), then A 1 = (I - D l M) x D [ and the 

proof follows as in the theorem. 
Let us now consider the error matrix A ' — X, and derive some bounds for its norm. 
THEOREM 1.2. Let X be a/z approximate inverse of A and let Y = I — AX. Then i/N(Y) < 1, 

N(A-i-X)-- SSSJ" 



Proof: 



so that 



l-N(Y) 

Y = I - AX 
AX = I-Y 

X- i A-' = (I-Y)- 1 

= I + Y + Y 2 + Y 3 + . . . 
A-*=X + XY(I + Y+Y 2 + Y 3 + . . .) 
A- , -X = XY(I-Y)- 1 

N(XY) 



N(A-'-X) 



l-JV(F)' 



where we have omitted the restriction that N(I) = 1 since it can be eliminated by the same process 
as used in the proof of Theorem 1.1. 

257 



COROLLARY. Under the same conditions as Theorem 1.2, 

N(X)N(Y) 



NCA-'-X) 



l-N(Y) 



PROOF: The proof follows immediately from the theorem and the properties of a norm. 

We have developed this latter upper bound as a corollary since, as we mentioned at the 
beginning of this section, this bound is frequently very much larger than the one in the theorem. 
It is obviously simpler and demands very little additional computer time but the information may 
be misleading (see the example at the end of this section). 

If we use the lefthand residual matrix, Y\=I—XA, the main difference from the above would 
be the need for N(Y\X) instead of N(XY). However, 

Y l X=(I-XA)X = X-XAX = X(I-AX)=XY 

so that actually the only difference comes from the use of N(Yi) instead of N(Y). Unfortunately, 
it is not true in general that Y\ = Y. In fact there are examples where there is a very large difference 
[4, p. 258]. So it is possible to have a good right inverse but a poor left inverse — good in the sense 
that the residual matrix is close to the zero matrix. Considering both residual matrices is of value, 
since they can give an indication of a poor inverse. 
Theorem 1.3. Let E = (I-XA)-(I-AX). Then 

NIA^-X)^^®-- 



2N(A) 
Proof: 

E={I-XA)-(I-AX) 

= (A-'-X)A-A{A-'-X) 

N(E) ^2N(A)N(A~i-X) 

that ME) 



2N(A) 

So the difference between the residual matrices gives us a lower bound which provides us with 
negative, but still useful, information. 

The question of which residual matrix to use in Theorem 1.2 and/or its corollary still remains. 
Using the theorem and requiring that the norm of both residuals be less than one, there should be 
very little difference. In general, however, the fact that the norm of one residual is small does not 
guarantee that the norm of the other will also be small: 

I-XA=A~ 1 (I-AX)A 

so that 

N(I-XA) ^N(A)N(A- l )N(I-AX) 
and similarly 

N(I-AX) ^N(A)N(A- 1 )N(I-XA). 

For an ill-conditioned system, N(A)N(A~ 1 ), the condition number, can be quite high. In practice, 
however, we find that the two residuals are almost invariably quite close, so that using Theorem 
1.2 with either residual should yield satisfactory results. Since inverses are calculated mostly in 
such a way as to almost guarantee the smallness oil — AX, it would seem that using/— XA in our 
error bounds would give more realistic information since we really expect a two-sided inverse. 
We also note that if it were particularly important to have I —XA small, it would be safer to invert 
A T , the transpose of A. 
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We have been considering only problems where the coefficient matrix had exact numerical 
entries. We would like to give now an error bound for the calculated inverse where there is some 
uncertainty in A [5]. Although we are working with the matrix^, we should really be working with 
A*, which differs from A by an amount satisfying N(A —A*) ^ a. We must, therefore, include this 
known difference in our bound for N( (A*) ~ ' — X) . 

Theorem 1.4. Let N(A - A*) ^ a and let R = 1 - A*X. Then i/N(R) < 1, 

N((A*)-> - X) < N(X)[N(I-AX)+aN(X)j 
N((A ) X)* 1 _ N(I _ AX) - a N(X) 

PROOF: We proceed exactly as in Theorem 1.2 and use R instead of Y, arriving at 

N(X)N(R) 



N((A*)-i-X) 



l-N(R) 



Now R = I - A*X = I - AX + AX - A*X = (I - AX) + (A - A*)X and N(R) ^N(I-AX) 
+ N(A — A*)N(X) *£ N(I — AX) + aN(X). Substituting this expression in the above inequality, 
we have the result stated in the theorem. 

Let us now turn our attention to the relative error of the difference between^ -1 andA". 

THEOREM 1.5. LetX be an approximate inverse of A and let Y = I —AX. 7/N(Y) < 1, 

1 - N(Y) 1 1 + N(Y) 



N(X) N(A->) N(X) 

PROOF: The proof follows immediately from Theorem 1.1. 

Corollary. i - N(Y) *= J^! X) ^ 1 + N(Y). 

THEOREM 1.6. Let X be an approximate inverse of A and let Y = 1 — AX. If N(Y) < 1, 

N(A-'-X) ^ N(XY) 1 + N(Y) 
NCA- 1 ) = N(X) " 1 -N(Y)' 

PROOF: The proof follows directly from Theorem 1.2 and Theorem 1.5. 

NIA- 1 - X) .._ 1 + N(Y) 

COROLLARY. N(A-) * N(Y) ' T^WO 

Actually, we have a simpler and better bound than this corollary. 

Theorem 1.7. Let X be an approximate inverse of A and let Y = I — AX or I — XA. Then 

N(A--X) * N (Y) 

Proof: A~ l — X = A~ l (I — AX) = (1 — XA)A~ l so that after taking norms, the inequaltiy 
follows. 
This bound is not necessarily smaller than the one given in Theorem 1.6. If N(Y) is small enough 
(say less than .1). the ratio in the above corollary can be expanded to the product 

[1 + /V(K)][1 +N(Y) +N(Y)t+ . . .] 
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or to 1 -f 2N(Y) if we ignore the higher powers of N(Y). Then the corollary becomes 



N(A-i-X) 
N{A-') 



N(Y)[1 + 2N(Y)] «7V(F), 



again ignoring the higher powers of N(Y). This indicates that in practice the bounds given in 
Theorem 1.7 and the corollary are quite close so that the bound given in Theorem 1.6 will probably 
still be the smaller one. 

Finally, we also have a lower bound. 

THEOREM 1.8. Let X be an approximate inverse, Y = I— AX, and let E = (I— XA) — (I — AX). 
Then i/N(Y)<l, 

N(A-*-X) ^ N(E) 1 - N(Y) 
N(A-i) ^2N(A)' N(X) 

PROOF: The proof follows immediately from Theorem 1.3 and Theorem 1.5. 

We conclude this section with a numerical example of the above bounds. Let N be the maxi- 
mum row sum norm and let A be the 6X6 segment of the infinite Hilbert matrix A, where 
aij = \j(i-\-j— 1). We chose this particular matrix because of its notoriously poor condition number, 
about 2.92 X 10 8 , but using only a 6 X 6 so that roundoff errors would not be too influential. Since 
we will use this matrix in all our examples, it might be good to give its exact inverse. For the rath 
order segment, if Ay t x = (typ), then 

(-l)' + J(n + i-l)\(n+j-l)\ 






l)[(i-l)!(7-l)!p(/i-0!(*-/)! 



At 



In our case (giving only the lower triangular part since it is symmetric) this turns out to be 

36 

-630 14700 

3360 -88200 564480 

-7560 211680 -1411200 3628800 

7560 -220500 1512000 -3969000 4410000 

-2772 83160 -582120 1552320 -1746360 698544 

We note that this is not quite the inverse of the matrix actually used in the examples. It is impossible 
to represent all elements ll(i+j— 1) exactly in the computer and hence the input matrix is slightly 
different. For further information on this point of inexact data of the Hilbert matrix, see reference 
12 (Computer Solution of Linear Algebraic Systems by Forsythe and Moler), section 19. We also 
indicate in section 6 how to avoid this problem. 

To calculate X we ran the program given in the above text (p. 68). The computer we used was 
a Univac 1108 and we used double precision for all matrix products. 

A word about notation. The standard computer notation .235 £ + 04 means .235 X 10 4 or 2350. 
The "1?" means exponent to the base 10 and the "+04" is the positive exponent. 
The absolute error: (Y = I — AX) 



(1) N(A-*-X) 

(2) N(A-'-X) 



N(X)N(Y) 
: 1-N(Y) 

. N(XY) 
1-/V(F) 
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= .197 £ + 06 



= .612 £-01 



(3) N(A-'-X) 



NJE) 

2N(A) 



.477 £-02 



The relative error: 



and for 



(1) /v(^') ** (y) i-/v(F)-- 168£ 01 



yv(^-'-Ar) 



#(7) 



= .162 £-01 



(2') ^^.^ •i±™=.522£-08 



N(A- 



N(X) 1-N(Y) 



N(A^-X) > N{El . 1-N(Y) = 393E _ Q9 
(6 ' N(A~>) "W{A) N(X) 66 ^ U 

N{X) -N(A->)- N{X) ' 



1+/V(F) ^" v " ' ~~l-N(Y) 
.117£ + 08 =£/¥(/*-') «.121£ + 08 



where A^" 1 ) = .11931731 £ + 08. 



Note: If we had been able to input the data exactly, the norm of the inverse would have been 
.11865420 £ + 08. 

It might be of interest to give the values of the previous quantities using the left hand residual, 
I-XA,foiY. 
The absolute error: 



(\)N(A"-X) 



N(X)N(Y) 
l-N(Y) 



.380 £ + 06 



(2) N(A-'-X) 



N(YX) 

l-N(Y) 



= .622 £-01 



(3) N(A-*-X) 



N(E) 

2N(A) 



= .477 £ -02 



The relative error: 



N(A-t-X) 1+N(Y) 

y ' N(A-*) ( ' \-N{Y) 



= .328 £-01 



N(A~ l -X) 
N(A~i) 



,N(Y) 



= .309 £-01 



N(A-'-X) ^ NjYX) 1 + N(Y) _ 
(2) N(A~<) ^InxT^NiY)' 



(3') 



N(A~'-X) N(E) 1-NjY) 



N(A->) 



2N(A) N(X) 
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= .387 £-09 



and for 



N(X) 



:N(A~ l ) 



N(X) 



1+N(Y) ^' v ~ ' ^1-N(Y) ' 
.116 £ + 08 ^A^- 1 ) ^ .123 £" + 08. 



We also note: 



7V(ZF) = .6024£-01 



N(YX) = .6024E-01 



and N(E) = . 2337 E -01. 



It will also be of interest to give the values of the previous quantities using a different computer, 
the CDC 6400, with Y=I — AX. This machine uses about twice as many digits in its operations 
as the Univac 1108 and this is why the results are so much better. We found that the values of the 
elements of X were integers (to 8 digits). 
The absolute error: 



(1) N(A- 


X) as v ' v ' 


= .323 £-01 


(2) N(A~ 


X) ^ N{XY) 


= .178 £-07 


(3) N(A~ 


x)>Tl 


= .555 £-09 



The relative error: 



N(A-i-X) < l+JTffl 



JV(/J-'-Z) 
Af(^-') 



;7V(F) 



.272 £-08 



= .272 £-08 






(3') *^^^ m^P = .468£-16 



N(A~ l ) 2N(A) N(X) 



and for 



N(X) 



■:N(A~ l ) 



N(X) 



1+N(Y) v ' l-iV(F) 
.1186542 £ + 08 ^iV^- 1 ) =£ .1186542 £ + 08. 



We also note that 



N(XY) = .178 £-07 and N(E) = .272 £-08. 
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2. Error Bounds for the Approximate Solution Vector of Ax b 

In this chapter we shall derive estimates for the closeness of the calculated solution vector 
to the true solution vector of Ax — b. Our approach will be to assume that we have available both 
the approximate inverse X and its residual Y = I — AX or I — XA. The necessity for this will* be 
explained by the quotations at the end of this section. 

Let z be the approximate solution vector of Ax = b and let r = b—Az, the residual vector. 

THEOREM 2.1. The bound for the norm of the error is given by 

N(A-'b-z>-™f, N(Y)<1. 

Proof: b—Az=r 

A- 1 b-z=A~ 1 r 

so that N(A-'b-z)^N(A-')N(r) 

. N(X)N(r) 

~~ i-yv(F) 

using our estimates for N(A~') from the previous section. 

Theorem 2.2. Let z=Xb. Then t/N(Y)< 1, 

N(XY)N(b) 

N(A <b-z) =s • 

1-/V(F) 

Proof: A~ l b-z = A~<b-Xb 

= (A~ 1 -X)b 

and N(A->b-z)^ N{XY)N(b) 

y ' 1-N(Y) 

again using information from the previous section. 

N(X)N(Y)N(b) 



Corollary. N(A~ l b — z) 



l-N(Y) 



We would like to make a remark stemming from the proof of Theorem 2.1. If the equality 
A~ ] r = A~ l b — z is ever used in practice, say for checking a process or computer program where 
the solution is known, we should note that we would be using X and not A' 1 . That is to say 

A-i r =A- 1 {b-Az)=A- 1 b-A- 1 Az=A- 1 b-Iz=A- 1 b-z 

but 

Xr = X(b-Az) =Xb-XAz = Xb-{l-Y x )z=Xb-z + Y x z, 

where Y, is the lefthand residual I — XA. Hence, we have an additional factor to keep in mind, 
Yiz. The need for this obviously depends on the quality of X; i.e., the closeness of XA to /. There 
are examples, however, where X may be a good inverse but where Xb may not be a good approxi- 
mation to the solution vector [1, p. 149]. So, in general, the correction factor Y x z can be of impor- 
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tance. To give some indication of just how important, let us give an example. 

In solving Ax= b where A was the 6 X 6 Hilbert segment and b was (l,2,3,4,5,6) r , we used Gaus- 
sian Elimination (with iterative improvement) to get X and the machine was an IBM 7094. For the 
maximum element norm, N(Xr) was of order 10~ 2 ; N{Yiz), of order 10 _1 ; N(Xb — z), of order 10 _1 . 
The need for the correction factor is clear. The fact that the Hilbert matrix is ill-conditioned results 
in the high order of the norms but the need is still there to get equality even for a well-conditioned 
matrix. Using a diagonally dominant tridiagonal matrix for A, we found that N(X r ) was of order 
10~ 7 ; N(Yiz), of order 10~ 8 ; N(Xb — z), of order 10" 8 . The need for the correction factor is again 
clear although not so vital. We would also like to emphasize that the residual to be used is the left 
hand one. As shown in section 1, there can be an appreciable difference in the two. 

The use of the lefthand residual (only) also gives us a bound for the relative error. 

Theorem 2.3. Let z = Xb. Then 

N(A-ib-z) 
N(A- J b) K } 

Proof: A- l b-z=A~ 1 b-Xb 

= (A~ 1 -X)b 

= (I-XA)A~ 1 b 

so that N(A- x b-z) Nil _v A \ 

N(A-'b) n ; * 

Once again this residual seems to be of more practical use than I— AX. (It was noted in the previous 
section that I — AX is ordinarily small by reason of the numerical procedure so that checking 
N{I — XA) is the important thing.) We can see that a small right residual/ — AX may not guarantee 
the smallness of the error vector from the following: 

A-'b-z = A-'b-Xb={A-'-X)b = A-'{I-AX)b. 

Hence, if A~ x has a large norm, the norm of the error vector may not be small. 

It can also be seen from the proof of Theorem 2.1 that a small residual vector may not indicate 
a good solution vector: 

A- 1 b-z = A~ 1 r 

so that N(A-ib-z)*zN(A-i)N(r). 

Again, if N(A~ l ) is large, N(A l b — z) could be large also. On the other hand, a large residual vector 
can indicate a poor solution vector. We show this in the following theorem. 
Theorem 2.4. 7/r = b — Az, then 

Proof: r=b—Az 

=A(A~ 1 b-z) 
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We would like to include here a brief discussion of an error bound for the computed solution 
vector of a system that does not have exact data [5, p. 151]. As in Theorem 1.4, let 

N(A -A*) *£= a, N(b-b*) ^ c, x*=A~% x**= (A*)'** and z=Xh, 

where X is the approximate inverse of^4. 

Theorem 2.5. Let R = I-A*X. Then i/N(R)< 1, 

N , _ **, N(X)[aN(z) + N(b-Az)+c] 
1 j ^ l-N(I-AX)-aN(X) 

Proof: z-x**=z- (A*)-^* 

= (A*)~*(A*z-b*) 

and N(z-x**) ^NKA^-^NWz-b*). 

Now A*z-b*=A*z-Az+Az-b + b-b* 

so that N(A*z-b*) ^N(A*-A)N(zY+N(Az-b)+N(b-b*). 

To get a bound for N[(A*)~ l ] we proceed as in section 1, starting with R =1 — A*X, which leads to 

in{ } J l-N(R) 

Also R = I-A*X 

= I-AX + AX-A*X 

and N(R) ^ N(I-AX) +N(A -A*)N(X). 

Combining all of the above, we have the result stated in the theorem. 

As mentioned at the beginning of the section, we have used the approximate inverse in all 
our bounds. This is most unfortunate because it means that we have to solve an additional n systems 
to get information about the solution of one system. It seems, however, that to expect any improve- 
ment in this situation is rather futile, according to some of the best men in the field. We shall 
quote from a few of them: 

A. M. Turing: "It is difficult to determine the accuracy of the solution of a set of equations without inverting the 

(1948) matrix." [6, p. 301] 

J. H. Wilkinson: "It does not seem possible to obtain a rigorous estimate of the accuracy of a solution to a set of equa- 

(1963) tions without some estimate for /V(y4 _1 ). . . . To achieve more . . . we need an approximate inverse 
and a bound for its error." [7, p. 126J 

E. L. Albasiny: "Unless one has some knowledge of the inverse of A we cannot immediately tell however how accurate 

(1964) our solution is." [8, vol. 1, p. 171 1 

L. Fox: "The search for the worst combination of uncertainties involves the computation of the inverse, and 

(1965) in solving linear equations we need a simpler and reasonably satisfactory criterion for the number of 
quotable figures in the calculated solution. If the coefficients vary widely in magnitude there is un- 
likely to be any satisfactory analysis . . ." [1, p. 159 J 
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A. S. Householder: "For error estimates it is advantageous to find at least an approximation to A 1 even when only a 

(1965) single vector x is required." [9, p. 91] 

W. Kahan: "The only disadvantage that can be occasioned by the lack of an estimate Z of A~ x is that there is 

(1966) no other way to get a rigorous errorbound for z — x." [10, p. 785] 

C. Moler: ". . . without a bound for the condition of A, no precise conclusion can be made regarding N(x n — x)." 

(1967) [11, p. 320] 

Hence, it would seem that the bounds given in this section are about the only ones possible. 
In some special cases, however, where A' 1 is easily obtained or estimated, some improvement 
in this situation might be hoped for. We will discuss some of these in section 5. 

3. Improvement Schemes for the Approximate Inverse 

In this section we will present three schemes for improving the inverse that we have found 
quite effective. Each has a different order of convergence but not without some expense. We then 
generalize to a scheme of any order of convergence. Next, we show how two iteration methods 
which are frequently used are special cases of our schemes. Finally we give a numerical example 
of the various bounds derived. 

The first scheme follows directly from the same procedure that yielded our error bounds for 
the approximate inverse in section 1. Let X be the approximate inverse and let Y=I—AX. Then 

A- l =X+XY+XY*+. . .. 

The righthand side of this equality gives us our first scheme: 

X k =X+XY+XY 2 + . . .+XY k , A^O. 

If N(Y) < 1, we have convergence as k goes to infinity and by taking enough terms we can get as 
close to the true inverse as we please. The scheme is ideal for computer execution, and double 
precision (at least) should be used for the products of the matrices and also for the sum. The 
result would actually be a double precision inverse and we have found this inverse to be very 
effective in yielding A~ l b, the actual solution vector for a linear algebraic system. 

The same scheme can give us error bounds for the closeness of X to A _1 . We have 

A~ 1 -X = XY+XY 2 +. . . 

so that we can get an excellent estimate for the difference just by taking enough terms. Hence, 
given any X we can tell how good it is, the only requirement being that N(Y) < 1. If that particular 
X is not good enough, we can take more terms from the series until we are satisfied: 

A~ 1 -X k =XY k + 1 +XY k + 2 + . . . 

= XY k ^(I + Y+Y* + . . .) 
so that 

N(XY k+1 ) 
MIA- 1 — **) ^-^^ — - 
In Ak) 1-N(Y) 

^ N(X)N(Y k+1 ) 

" i-yv(F) 

N(X)N(Y) k+1 



l-N(Y) 



As mentioned previously, the first bound is by far the best. However, the last bound can be quite 
useful. Having calculated X and F, we can easily determine k to satisfy our requirements on 
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N(A~ X —X) (realizing that this is definitely an overestimate), and proceed to sum the correct number 
of terms from our series and also determine the kind of multiprecision that may be needed to do the 
calculation correctly. Using this last bound also eliminates one matrix multiplication, to getXY k+1 . 
If we use the lefthand residual, the procedure is essentially the same and we have for 

N(I-XA)=N(Y)<1, 

A-'=X+YX+Y 1 X+. . . 

X k =X+YX+Y 2 X+. . . + Y k X 

and 

N(Y k+l X) 

_ N(X)N(Y k+l ) 

* l-N(Y) 

^ N(X)N(Y) k+ > 

* \-N(Y) 

The second scheme we are interested in is very similar but the residual matrix must be recal- 
culated at each step. We start as before but now let X a =X and define a recursion formula: 

X k =X k - 1 +XY k =X+XY+XY i + . . .+XY" 
Y k = I-AX k , k>0. 

If N(Y) < 1, we are assured of convergence as k goes to infinity and with Y k as defined we proceed 
exactly as in section 1 to get 

in our derivation for this bound we assumed X^ x exists. It is clear that one of our basic assump- 
tions from the beginning has been that X~ x exists. It seems, however, that we should be more 
explicit about the existence of X~\ 

X k = X+XY+XY* + . . .+XY k 
=X(I+Y+Y 2 + . . , + Y k ). 



Now 
Hence 



(I-Y)-i = I + Y+Y 2 + . . . + Y k + Y k+1 + . . . . 

I + Y+Y- + . . . + Y k =(I-Y)~ l -Y k+l (I + Y+Y i + . . .) 
= (I-Y)-i-Y k+l (I-Y)- i 
= (I-Y k + 1 )(I-Y)-K 

Since N(Y k+ ') ^N(Y) k+l ^N(Y) < 1, (I—Y k > ')-' exists and hence we are assured of the existence 

ofZ,7'. 
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In this scheme the residual actually gets multiplied by itself as the process continues: 

Y i = I-AX l 

= I-A(X+XY) 
=Y-AXY 

= (I-AX)Y 
= F 2 . 



Y k =I-AX k 

= I-A(X k - 1 +XY k ) 

= Y k . 1 -AXY k 

= Y k -AXY k (by induction) 

= Y k+l 

Since N(Y k ) ^ N(Y) k+t < N(Y) , our condition for convergence and the validity of the error bound 
is still N(Y) < 1. 



Now 



X k =X+XY+XY 2 + . . .+XY" 



so that 
Since 

we have 



l-p A 



1-P 



N(X k ) ^N(X)[l + N(Y)+N(YY- 



-=l+p + p 2 + . . .+p A ', p^l, 



N(X k )*zN(X) 



\-N(Y)'^ 
1-N(Y) ' 



■N(Y) k ]. 



We also have 



X k Y k =(X+XY + XY 2 +. . .+XY k )Y k +* 



:XY k+1 + XY k+2 + . . .+XY 2k+l 



so that 



N(X k Y k )^N(XY k+1 )[l+N(Y) + . . . + N(Y) k ] 



--N(XY k+l ) - — tllLl — 



Hence 



N(X k Y k ) < N(X k Y k ) 
l-N(Y k ) ^1-AW +1 



^. N(XY k+1 ) 

* 1-/V(F) 
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and the best bound for our second scheme is better than that for our first. 
Summarizing the bounds for the second scheme, we have for N(Y) < 1, 

N(A ' X.) - N(XkYk) 

~~ l-N(Y)'- + i 
^ N(X)N(Y) k+1 

* i-yv(F) • 

It will be noticed that the last bound is exactly the same as for the first scheme, so that any improve- 
ment must come by using the first bound. 

In spite of the fact that we recalculate the residual matrix at each step, we used it only in the 
error bound expression. It would seem useful to try to get it involved in the calculation of the new 
approximate inverse also. We can do this quite simply and very effectively by returning to the basic 
relation at each step. The process would be as follows: 

A~ i =X+XY+XY 2 + . . . 

with convergence if N(Y) < 1. 

X X =X + XY=X(1+Y) 

Y { = I-AX X = Y 2 . 

Using this expression for Y u we repeat: 

A~ i =X l +X 1 Y 1 +XiY? + . . . 

with convergence if N(Yi) < 1 
but N^) ** N(Y) 2 . 

X 2 =X 1 +XiY 1 =X i (I+Y 1 ) 

Y 2 = I-AX 2 = Yj = Y 4 = Y 2 \ 

Xr+l=Xr(I+Y r ) 
I r+i = ' AX r + 1 . 

Y r+ i=I-A(Xr(I + Y r )) 

^ I r /I A. yY r 

= Yi 



In general we have 



Now 



that 



Y r = Y* r . 



Hence, if N(Y) < 1 we have convergence for our third scheme and since the error matrix is squared 
at each step, we have quadratic convergence. 
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It is interesting to observe that this third scheme actually amounts to taking a very large num- 
ber of terms in our first series at each step. For example, 

X i =X 2 (I+Y 2 ) 

=X 1 (I+Y 1 )(I + Yf) 
=Z(/ + F)(/ + F 2 )(/ + F 22 ) 
=X+XY+XY* + . . .+XY 7 



and in general, 



* r =X(/+F)(/+F 2 ) . . . (Z + F 2 '- 1 ) 
=X+XY+XY* + . . ,+^F 2 ''- 1 . 



The recurrence formula, however, is much simpler and more efficient. 

The bounds for the third scheme follow directly as before. If N(Y) < 1, 

NiA-'-X,.)-"^ 



l-N(Y r ) 

< N(X r )N(Y r ) 

" l-N{Y r ) 

< N(X r )N(Y)* r 

l-N(Y) 2r 

^ N(X)N{YY r 
* \-N{Y) 

The existence of X? 1 can be shown in the same way as in our second scheme. 

If we had used the lefthand residual matrix, the obvious changes would follow as before and 
we also have Y r X r =X r Y r . 

We can generalize our approach in the sense that we can define a recursion such that the resid- 
ual (error) matrix can be any power of the previous residual matrix. For example, to get third- 
order convergence, we take the first three terms of our basic series and repeat after each step. Thus 

X s+1 = X s (I+Y s + Y-i) 
i«+i — I AX s+ \ 

with Xq = X, Yq=Y=I — AX and N(Y) < 1. The proofs follow exactly as in the second-order case 
and we have 

y —yzs+i 

However, another matrix multiplication needed at each step and the possible loss of significant 
figures are factors to be considered in the practical application of this scheme. 
For any order convergence, say t, we have 

X s+l =X s (I+Y s +Y* + . . . + FJ" 1 ) 
Fs+i = / AXs+i. 

The proof is straightforward by induction on t and we also have 

F S+1 = F' S+1 . 
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The error bounds follow directly as before: 
If N(Y) < 1, 

N(X,Y.) 



N(A '-*,) 



1-N(Y,) 

^ N(X,)N(Y.) 
~~ l-N(Y,) 

\-N(Yy x 
i-yv(F) 

the first being the best, but the last being useful in determining s as was indicated on page 266. 

We note that the multiprecision needed for the effective use of these higher order convergence 
schemes might be a minimal problem in view of the hardware being developed by the computer 
industry. At the end of section 1 we gave an example using different computers and we found that 
the CDC 6400 gave extraordinary results. 

Let us now discuss the relative error of the approximate inverse derivable from these schemes. 
For each scheme we can use 

A-i-X k =A-HI-AX k )=(I-X k A)A-* 

so that 

N(A~<) y A '- 

where Yi ; may represent the lefthand or the right hand residual matrix. Hence, lor our first scheme: 

N(A-i-X k ) A;/l/ , 









N(A-*) ^' V(J '- 


for the second: 






N(A->-X k ) 

N(A-*) n ' 
=£AT(y)*+i; 


and for the third: 






N(A-i) n ; 
^N{Y)* r . 


In general for tth 


order 


convergence: 








MA- 1 ) K ' 

^N(Y)' S . 



Hence, if we use this relative error as a criterion for stopping the process, it is evident that the 
third is far better than the first two. There are also other considerations which will be mentioned 
shortly. To get a better bound we can use the relation 

I 1 + /V(F,,) 



N(A~<) N(X k ) ' 
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which follows exactly as in section 1. Hence, for our first scheme: 

N{A-!-X k ) ^ N(XY*+*) 1 + N(Y) 
N(A-i) '"■ N(X) ' l-N(Y) 

V ' \-N{YY 



for the second scheme: 



N(A-*-X k ) ^_ N(X,.Y k ) l + N(Y k ) 
N(A-i) ~~ N(X k ) l-N(Y k ) 



N(Y k ) 



l + N(Y k ) 
l-N(Y k ) 



1 + ]\f(Y)k-+i 
n ' l-JV(F)*+i' 



and in general for the rthe order convergence scheme: 



N(A~ l - 


-X.) 


„ N(X.Y.) 


N(A~ 


') 


' N(X S ) 
^ N(Y.) ■ - 



■:N(Y) fS • 



1+N(Y 8 ) 

l-N(Y s ) 

1+N(Y 8 ) 

l~N(Y s ) 

i + N(Yy 8 
i-N(Yy 8 ' 



Let us compare the three main schemes we have discussed. Our basis for comparison will be 
the number of new matrix multiplications needed to advance one step in the process and get the 
error bound for that particular iterate, presuming the previous one was not small enough. This 
means, of course, that the error bound must be calculated at each step. We have written the third 
scheme as X k +X k Y k instead of X k (I+Y k ) since the product X k Y k was calculated in the previous 
step for the evaluation of the error bound and is available. We will identify the first, second, and 
third schemes by Scheme I, Scheme II, and Scheme III in the following table and list the quantities 
that must be calculated in going from one step to the next. 





Scheme I 


Scheme II 


Scheme III 




Step k Step k+\ 


Step k Step k+\ 


Step A- Step k+1 


Iterate 


Y Y 

N{XY k+ i) N{XY k ^) 
l-JV(F) l-N(Y) 

1 
XYk+2 


X k X k +XY«+* 

Y k Y k Y 

N(X k Y k ) N(X k+1 Y k+1 ) 
l-N(Y k ) l-N(Y k+1 ) 

3 

XY«+\Y k Y,X k+1 Y k+i 




Residual 


Y k Y k Y k 


Bound for N(A~ l -X k )... 
Additional products 


N(X k Y k ) N(X k+1 Y k+l ) 
l-N(Y k ) l-N(Y k+l ) 

2 

Y k Y k , X k+i Y k+l 



We see that the second scheme has one more matrix multiplication than the third and, since 
the latter has a much better rate of convergence, the third scheme is definitely better than the 
second. In comparing the first and the third, the question of the additional matrix multiplication 
to get better convergence arises and cannot be answered in general. The exponent of Y in Scheme 
I may be such as to demand multiprecision hardware and this may be worse than performing an 
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extra matrix multiplication. In some cases too, this extra multiplication may itself demand multi- 
precision hardware. In view of the quadratic convergence, though, many fewer steps should be 
required for the third scheme and probably multiprecision hardware would not be needed (beyond 
double precision, that is). We note that in our second scheme, the convergence is better than in the 
first and we have proven that the bounds are also better. Two additional matrix multiplications 
may not be worth it, however. Considering all of these factors and also keeping in mind the relative 
error, we feel that the third scheme is the most practical. 

We would like to show now that two frequently used iteration procedures to find the inverse 
of a matrix are actually included in our development. The first is defined by the recurrence 

X k+1 =X+(I-XA)X k9 *2*0, (A) 

with Xo = 0. Letting Y/ X = I — AX^, we can rewrite this recurrence as 

Xk + 1 — X ■+■ Xic — XAXk 

=X k +X(I-AX k ) 

=X k +XY k . 
Now 

Y k =I-AX k 

= I-A(X k . 1 +XY k - 1 ) 

= I-AX h - ] -AXY k - l 

= (I-AX)Y, , 
= /T,,-, 

Y k =Y k . 

X k+1 =X k +XY k , k^Q 

Yk+] = I ~ AXk + \, 

which is exactly our second scheme (the starting value is the only difference). However, in the origi- 
nal form, only the lefthand residual can be used whereas in our development either one is appli- 
cable. We note also that for the convergence of (A) each eigenvalue of the righthand residual must 
be of modulus less than one. Although this seems to imply the necessity of working with both 
residuals, we realize that this is not the case since the two residuals are similar: 

I-AX=A{I-XA)A~ l . 

The second frequently used iteration is defined by 

X k+1 =X k (2I-AX k ), k^O (B) 

where X is arbitrary. The immediate difficulty is the choice of X and hence this iteration is more 
useful as an improvement scheme. If we let X = X and Y k = I — AXk, we can rewrite (B) as 

X k+1 =X k (I+Y k ) 

which is identical to our third scheme. The changes to make use of the lefthand residual, once (B) 
is written this way, follow exactly as in our development. 
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and by induction we have 

So this recurrence can be written 



So we see that both of these iteration formulas can be thought of as having their foundation 
in our basic relation that 

A-'=X + XY+XY*+ .... 

We have repeatedly mentioned that in our schemes we can improve either the lefthand or 
the righthand inverse. This is a very important consideration in the solution of Ax=b, as will be 
pointed out in the next section. 

For a numerical example, let A be the 6X6 segment of the Hilbert Matrix and to generate 
the approximate inverse, X, we use the same program and computer as in section 1, except that 
we used no iterative improvement scheme. We used double precision accumulation for the products 
XY and our criterion for stopping was when the elements of XY k were in magnitude less than 10 ~ 8 . 
The final results for each scheme were identical. 

For convenience we list the results separately and then make some comments about them. 



Scheme I 

N(A- 1 -X)*zN(XY+XY 2 + . . . +Xy 5 ) =.412 £ + 05 

N(A-*-X 5 )^f^^ =.227 E- 01 

^ N(X)N(YY 
" 1-N(Y) 

N(A-i-X) 



= .199 £-02 



yv^- 1 ) 



N(Y) =.234 £-01 



and for 



N N(A-* ) ^ N(XY+XY2+ ■ ■ ■ +ZF5) *^zP =-360 £-02 

N{A-*-X s ) ^ N(XY°) l + N(Y) = 
N(A-i) ~~ N(X) l-N(Y) 

^ N W 6 \-N ( (Y) =.170 £-09 



l + N(Y) v ; l-yV(F)' 
.117 E + OS^NiA- 1 ) s; .122 £ + 08 
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and for 



and 



Scheme II 
N(A-^-X,)< 1 2 ' " = .114£-04 



1-N(Yk) 

^ N(X)N(Y)« 
' 1-N(Y)« 



= .502 £-04 
= .194 £-02 



N(A-t) ~ N{Yi) =.421 £-11 



NjA-i-Xt) „ N(X S Y 5 ) 1+JV( Y,) 9M E _ i2 



/V(/M) ' tf(Z 5 ) i-yv(F 5 ) 



l + N(Y) e 



NiXi) -,n<4-*)~ N{Xi) 



i+N(Yt) v ' i-yv(n)* 

.11931731 £ + 08=syV(/f-i) «. 11931731 £ + 08 

Scheme III 
N(A-i-X a )< Y&£L =.840 £-05 



= .505 £-04 
= .108 £-05 



^ N(X 3 )N(Y 3 ) 

" i-iv(y 3 ) 

< /V(X)/v(Fr 

1-/V(F) 

*^£f«W.> = .424£-H 

#(/*-' -X 3 ) ^(* 3 F 3 ) 1+7V(F 3 ) 7ruF -i9 
AT(^-i) ^ /V(Z 3 ) 1-JV(F 3 ) 

1 + /V1F) 8 
*W>''IZ]5<F)i =.886 £-13 



^^ .JV^uJ^ 



1 + /V(F 3 ) " v * ; 1-/V(F 3 )' 
.11931731 E + 08^N(A~ l ) *£ .11931731 £ + 08. 
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We would like to make some comments about this example. One of our theoretical results 
was that 

N(Y k ) **N(Y) k +K 

NowTV(F) 6 ~ .162 £ — 09, so that this inequality does actually hold. 
Another inequality we proved was that for Scheme II 

N(X k Y k ) ^N(XY^) 



l-N(Y k ) l-N(Y) ' 

We notice here, however, that this is not the case for k = 5. In investigating further we found that 
this inequality did hold for k= 1, 2, 3, and for A : = 4 there was almost equality. In examining the 
significant digits involved, we found that after A = 3, at which point N(Y 3 ) ~-l £ — 08, we had 
gone beyond the accuracy of the numbers with which we were dealing. We were working in double 
precision with numbers originally of the order £ + 08 and have now gone beyond 16 digits. Hence, 
any information beyond this point becomes meaningless. 

The difficulty is more apparent in dealing wit\\X k Y k andZF fc+1 than with just Y k and Y because 
of the matrix multiplications involved which lead to more serious round off errors. This is the 
reason that we mentioned in this chapter the possible need for multiprecision when working 
with these schemes. 

The contradiction in Scheme III is different but comes from the same source. Here 

N (X 3 Y 3 ) > N(X)N(Y)* 



l-N(Y 3 ) l-N(Y) 

and A = 3, the quadratic convergence making the difficulty show up earlier. 

This example very clearly illustrates the point of the question asked in the Introduction: 
How do we know these numbers are the solution of our problem? The bounds developed here 
certainly help in answering this question. 

4. Improvement Schemes for the Approximate Solution Vector of Ax= b 

Let us now examine ways of improving the solution vector, z, of a linear algebraic system, 
Az—b. As in the previous section, the residual, r—b—Az^ plays an important role. The standard 
approach is to solve the system Ax=r whose solution, d, is an improvement to be added to z. 
The usual convergence proof depends on Wilkinson's backward analysis which says that z is the 
exact solution of a slightly perturbed system (A + E)x = b. The condition for convergence of the 
iterative improvement scheme [12] is that N(A~ l E) < i. We would like to present a slightly different 
scheme and give two convergence proofs for it. We will then give a different convergence proof 
for the standard approach, one we find more meaningful. 

The alternate scheme differs in that the residual that is used comes from the system just 
solved and not from the original one; that is, r =b—Az and if do is the solution of Ax = r , then 
ri= ro — Ado. In general we have 

r k +i=r k —Ad k 

and dk is the approximate solution of Ax = r k . 

THEOREM 4.1. Let the vectors r k and d k be defined as above. Then i/N(r k + i) < aN(r k ),0 < a < 1, 
r k approaches zero as k goes to infinity and z+ do+ di + . . . + d k converges 
to A-!b. 
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PROOF: We have 



r {) = b — Az 
r x = r () — Ado 
r> = fi — Ad\ 



that 



or 



Nc 



rjfe+i = r k — Ad k 

r +ri + . . . + rjk + i=6 + r +ri+. . . + r A --/4(z+do+di + . • . + <4) 

rfc+i=6— A(z + do-+di+. . , + d k ). 

N(r k+1 ) < aN(r k ) <aW(r,_,) < . . . < a*+W(r ). 



Hence, as A: increases to infinity, a k+1 goes to zero, and then N(r k +i) goes to zero which implies 
that />+i goes to zero. Therefore A(z + do+di + . . . + <4) approaches 6 which means that z+rfo 
-\-d\-\- . . . +A- converges to A~ x b. 

Instead of summing the residuals, if we use Wilkinson's backward approach, we could say 
that d k is the exact solution of a slightly perturbed system and we would have 

(A+E k )d k =r k 

where we have subscripted the E to show the dependence on each system. With the introduction of 
E k we have another convergence criterion. 

THEOREM 4.2. Let the vectors r k , d k and the matrix E k be defined as above. Assume N(E k ) < c 
for all k. Then if N(d k+/ ) < pN(d k ), < p < 1, r k approaches zero as k goes to 
infinity and z + d + di + . . . + d k converges to A _1 b. 



Proof: 



but 

so that 
and 
Now 
Hence, 



(A-\-E k )d k =r k 

Ad k + E k <l k =r k 

E k d k = r k — Ad k 

r k+ i= r k — Ad k 

r k+ i = E k d k 

N(r k+l )^N(E k )N(d k ). 

N(d k )<pN(d k - 1 )<p*N(d k - 2 )<. . .<pW(rf ). 

N(rt+i)<cp*N(do) 



so that as k goes to infinity, p k approaches zero and r k +i approaches zero and the conclusion 
follows just as in Theorem 4.1. 
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COROLLARY. Under the same conditions as the theorem, we have N[A _1 b — (z+ do+ . . . 



Proof: 



+ d k )]^^N(d ). 



A- 1 b—(z + d +. . ,+d k ) =dk+i + dk +2 +. . • 
N[A-*b-{z+do+. . . + d k )]^N(d M )+N(d k +*) + . . . 

^p*+W(do)+P* +2 N(do) + . . 

= P*+ 1 (1+P + P 2 +- . .)7V(d ) 
D fc+1 

-iV(A). 



1-p 

We will see that these two conditions are actually contained in our next theorem. There 
N{I — XA) < 1 is our condition, and this is the p here; and also N(I — AX) < 1 is used, and this is 
the a here. 

It would seem, however, that repeatedly calculating the residuals from the primary equation, 
b—Ax, is a better scheme since the magnitude of the residuals that are being used may be much 
smaller for the alternate method. Using 'the smaller numbers might demand multiprecision rou- 
tines, more than the double precision necessary for the calculation of each residual in the standard 
approach. On the other hand, it has been found that N(rjc+i) is not always smaller than N(rjc) for 
the standard scheme, and in similar problems using the other scheme we have found that the resid- 
uals do decrease at each step. This could be due to the fact that in the second scheme we make one 
addition at the end, z + do+di + . . . -\~dk, for our final answer and thus only one rounding error 
occurs. In the standard scheme, however, this is done at each step and if this rounding error is 
included in the convergence proof, we would find that the error, ^ _1 6 — #A-,does not approach zero 
but only gets small [12, p. 110]. Hence, this convergence proof is really for what Forsythe calls an 
"Almost Theorem" [13, p. 9], especially since Kahan has found a counterexample to the standard 
approach [10, p. 78]. In practice we have found both schemes produce identical results. 

Let us now turn our attention to a different convergence proof for the standard scheme. It 
will be useful to write out a few steps of the process. Let du be the computed solution at each step. 



(1) solve Ax = b: 



(2) solve Ax = r\ 



z = Xb 

r=b — Az 
= b-A(Xb) 
= {I-AX)b 

d =Xr 

= X(b-Az) 
= Xb-XAz 

= (I-XA)z 

= Yz, Y=l-XA. 
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improved solution: 



and 



(3) solve Ax= n 



improved solution: 



In general, we have 



z\ = z -f do 
= z+Yz 
= {I+Y)z 

/i = b — Az\ 

= b-A(z + d ) 
= b—Az—Ado 
= r-AXr 
= (I-AX)r. 

di=Xri 

= X(I-AX)r 
= (I-XA)Xr 
= YXr 

= Y 2 z 

z-i = Z\ + d\ 
= (I+Y)z+Y-z 
= (I+Y+Y*)z 

n =l> — Az 2 

= b-A(zi + di) 
= b—Azi—Adi 

= r 1 -AXr 1 

= (I-AX)n 

= (I-AX) 2 r. 

dk-i = Xrk-i 
= Yd k - 2 

= Y«z, 

Zk= Zk~i + dk-i 
= (/+K+F 2 +. . . + Y*)z, 

r k =(I-AX)r k -i 

= (l-AX) k r 
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406-680 0-71-4 



for by induction 

=X(I-AX)*r 

=X{l-AX){I-AX)k~h 
= YXr k -i 
= Yd k -x 

= F*+ 1 z; 

Zk+1 = Zk + d/c 

= (/+F+F 2 +. . . + Y*)z+Y k+ h 
= (/+F+F 2 + . . . + F*+F*+»)z; 

rk+i = b — Azk+i 
= b-A(z k +d k ) 
= b—Azk — Adk 
= r k — AXr k 

= (I-AX)r k 

= {I-AX)«+ i r. 

Theorem 4.3. Let the vectors d k , z k , r k and the matrix Y be defined as above. Then i/N(Y) < 1, 
z k converges to A _1 b as k goes to infinity. 

Proof: z k =(I+Y+Y* + . . . +Y k )z and, as we proved in Lemma 1.2, if N(Y) < 1, /+F+F 2 + 
. . . +Y k converges to (/— F) _1 as k goes to infinity. Also 

Y=I-XA 

XA = I-Y 

A- l X-*={I-Y)-K 

Hence, as k goes to infinity, z k converges to 

(I-Y)- 1 z=(A- l X- 1 )(Xb)=A-*b. 

Corollary 1. Let d k and r k be defined as above. Then if N(Y) < 1, d k goes to zero as k goes 
to infinity and i/N(I — AX) < 1, r k also goes to zero. 

Proof: 

dk-x = YH 

so that 

N(d k .i)^N(Y)"N(z). 
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Since N(Y) < 1, the righthand side goes to zero as k goes to infinity so that the lefthand side does 
also, which implies that dk-i goes to zero. Since 

r k =(I-AX)*r 
and 

N(r k ) ^N(I-AX) k N(r), 

if N (I — AX) < 1, the same argument implies that r k goes to zero as k goes to infinity. 

We still have linear convergence but we feel using N(Y) < 1 is more meaningful than 
N{A~ 1 E) <i. Since we need X if we are going to use any bounds for N(A~ l b — z k ) , it is also a 
practical criterion. 

From this theorem we can also get an estimate of the improvement at each step as well as a 
final error bound. 

COROLLARY 2. With the same definitions and conditions as in the theorem, 

(a) N(z k - z k _ t ) *£ N( Y) k N(z) and 



Proof: 



that 



so that 



N(Y) k+1 
(b) N(A-*b-z k )^y^^N(z) 



(a) z k — Zk-i = d k -\ 

= Y*z 

N(z k -z k - l )^N(Y)m(z). 

(b) A-ib-z k =(I-Y)-h-z k 

= (/-F)- 1 z-(/+F+ . . . +F*) 
= (Y k ^+Y«+*+ . . .)z 

N(Y k+l ) 
N(A-ib-z k )^ f^^ N(z) 

N(Y) k+1 
l-N(Y) nz) ' 



This is a smaller number generally than that bounding the difference between two successive 
steps in the iteration. In fact, if N(Y) ^ i, 



l-N(Y) 
and 

N(Y) k+1 

^N(Y)m(z). 

We note that for the convergence of the sequence of zk we need only that the norm of the left- 
hand residual be less than one. This implies that the residuals, b—Azk, need not go to zero as k 
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increases and this has been the case in some problems as previously mentioned. This fact also 
emphasizes that in the solution of Ax = b, the lefthand inverse is the more important one. As was 
pointed out before, a good right inverse may not imply that the error in the solution vector is small. 
Hence, if we intend to use an approximate inverse in our solution of Ax = b it is very important 
that the improvement scheme used be one that improves the lefthand inverse. This is very easily 
done by the methods of the previous section. 

5. Error Bounds for the Approximate Inverse of Special Matrices 

In this section we will briefly discuss some results concerning a symmetric matrix, derive an 
error bound for the inverse of a triangular matrix, and then investigate the LU decomposition of a 
matrix A for possible error bounds on^ _1 which we will apply to the solution of Ax = b. 

Since we will be developing results that apply to some special matrices, it might be good to 
give some properties of these matrices so that they will be readily available. We will not give any 
proofs as these can be found in various texts on linear algebra or numerical linear algebra as listed 
in the bibliography. 

Let A be a real symmetric matrix. 

1. The eigenvalues of A 2 are the squares of the eigenvalues of A. 

2. The eigenvalue of maximum modulus of A*A is the square of that of A (y4* is the con- 
jugate transpose of A and hence, in our case, is just the transpose of A, A T .) 

3. The spectral norm of any matrix A, N g (A) , is defined as the square root of the eigenvalue 
of maximum modulus of A* A. Hence, in our case, N S (A) is just the absolute value of the eigenvalue 
of A which has maximum modulus. 

4. A is orthogonally similar to a diagonal matrix; i.e., A = T DO where is orthogonal, 
T =0~ 1 and D is diagonal. 

We will also use a few general relations which we will mention at this time. 

1. Let N/(A) be the Frobenius norm of A. Then if A and B are orthogonally similar, N/(A) 
= Nf(B). (See section 1, where we proved this for unitarily similar matrices.) 

2. We recall a theorem of Schur (also used in sect. 1): If A is an arbitrary matrix, then it 
can be transformed to a triangular matrix by means of a unitary matrix; i.e., a unitary matrix U 
exists such that U*AU= T, where U* = U~ 1 and T is upper triangular. 

We will now define a relation between matrices which may not be too familiar. 

DEFINITION. A is dominated by B if and only if | a^l ^ by for all i and j. Symbolically, we write 

A <^ B. Some basic properties of this relation are that if Ai <^ Bi and A 2 <^ B2, 

then Ai + A 2 <^ Bi -f B 2 and Ai A 2 <^ BiB 2 . 
Finally, we will use the symbol/ to represent the n X n matrix whose elements are all ones, and 
U will be used for the eigenvalues of the matrix under consideration. 
With this background we can proceed to our first theorem. 

THEOREM 5.1. If A is a real symmetric matrix and is O < r ^ |1 4 |, then N S (A _1 ) ^ 1/r. 
Proof: r^ \l { \ 







1 


1 
^ — 

r 




\li(A~ 


l )\ 


^1/r 


max 


\h(A~ 


l )\ 


^1/r 




N S (A 


- 1 ) 


sU/r. 



In our next few theorems let the following notation hold. Let X be the approximate inverse of 
A and let R = OXO T = (r#) where is an orthogonal matrix such that A = T D0. 
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THEOREM 5.2. Let A be a real symmetric matrix. Then if 



Proof: 



and 



Now 

so that 
Hence 

Now 



N/(I-AX)«e, N,(A-'-X) ^ n 1 ' 2 

A- 1 -X=0 T D~ 1 0-0 T RO 
= T (D~ 1 -R)0 

N f (A- 1 -X)=N(D- 1 -R) 

= N[D~ l (I-DR)]. 

I-AX = I-0 T DRO 
= 0'(I-DR)0 

Nf(I-AX)=N f (I-DR). 

N f (A~ l -X) *N f (D-i)N f (I-AX). 



min |lj 



so that 

and if N/(I—AX) ^e, we have 



D'< 



N f (D-*) 



1 



min |//| 



/ 



min 



N f (A^-X)^n ] 



min I //I 



COROLLARY 1. Under the same conditions as the theorem except that if I — AX <^ kj, we have 

N y -(A-i-X) ^n' j / 2 . k , 1 , • 
min |li| 

PROOF: Since N/(J) =n, we use nk instead of e as in the theorem. 

Corollary 2. Let all the eigenvalues of A be the same, say a. Then 



N f (A^-X)^~ 
\a\ 

Proof: We use the fact that D= al and the property of a norm. 

Theorem 5.3. Let the n eigenvalues of A be distinct. Let E= (I — XA) — (I — AX) and let 
E, = OEO T =( eij ). ThenifA=\^-x tt , 



N/(A-'-X)^n max(|di|, 
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L-li 



Proof: 



E=(I-XA)-(I-AX) 
=AX-XA 
= T {DR-RD)0. 



If 



E x = OEO T , E 1 = DR-RD. 
In terms of elements, this last equation can be written 

= , i=j 



or 



lJ li-lj 



, i^j 



and for i=j\ ry is not determined but they should be close to l\ 1 . Now consider D~ l —R. The diago- 
nal elements are d\ and the off-diagonal, — e^\{l\ — lj). Then 



(«• Hi 



so that 



D~ l -R< max |A|, 



N/iA-i-X) =N f (D- 1 -*R) 



n max I |d*|, 



The spread of the eigenvalues of a matrix can give a good idea of the condition of a matrix; in fact, 
the ratio of the modulus of the largest to that of the smallest is defined as the P-condition of a matrix. 
We see from our error bound that a large spread does not necessarily mean we cannot get a good 
approximate inverse. This bound implies, however, that if two eigenvalues are very close, we 
might expect trouble but that this might be counteracted if E is small enough. Hence, once again 
knowledge of the lefthand residual matrix is useful. 

Corollary. /FE = 0, NfiA- 1 - X) ^ n max | d 4 1 . 

i 

PROOF: If E = , E x = and then nj = for i # /. Hence 

D 1 -R = diagonal (di) 
< max | di | J 

NfiA- 1 —!) ^n max \d t \. 



that 



Let us now consider another special case, a triangular matrix T, and derive an error bound 
for its inverse. Let N rs be the maximum row sum norm. 

hi 



Theorem 5.4. Let T be a triangular matrix and let e = max 

i 

(l + e)"- 1 



tu 



for all j. Then 



N rs (T- 



min | tu | 
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Proof: Let T=D-U=D(I-D- l U)=D(I-E), D diagonal, U strictly upper and 
E = D- i U=(e iJ ) where e«=fy/f«. Then T~>= (I-E)-'D- 1 . Let us examine (/-£)-». Since £ 
has zero diagonal, £'< = and (/-.E')- , = / + £ + £ 2 + . . . +£»-». Let e = max |e y |, for all;. 
Then we can write E < eP where ' 



P = 



1 




Let us examine the various powers of P. 



pi = 






1 


2 


3 


4 . . 


/i-2 








1 


2 


3 . . 


n-3 



2 
1 





P* = 



1 3 6 10 15 



71-2 

2 



10 
6 
3 
1 
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In P 3 we have written 



71-2 



instead of 1 -h 2 -h 3 + . . . +n — 3 which comes from multiplying 



the first row of P times the nth column of P 2 . We also note that n — 2 can be written 



("T 2 ) 



fhich 



we will use in P 2 ; of course I n 1 = 1, the element in the (1, n) position of P 1 . Generalizing this 
pattern and realizing that P n =Q, we find thatP n_1 has all zeros except for a 1 in the (1, n) position 



which we will write as 



n-2 

n-2 



Let us write the first row (only) of the next few powers: 



Pi=[o 











1 


4 


10 


20 


■ ("i 2 )] 


p?=[o 














1 


5 


15 


■• (7 2 )] 



We restrict ourselves to the maximum row sum norm and our interest is in 

N(I + E + E 2 + . . . +£>-!). 

As mentioned in our first chapter, it will be better to add the matrices first and then take the norm. 
It is clear that the maximum row sum will come from the first row. Let us rewrite the first rows of 

these matrices, then, using the binomial expression 



p\- 


= 


(1) (1) (I) C 


p\- 


.[. 


» (!) ( 2 ) (! 


p\- 


.[. 


• • ( 2 ) (I 


p\ 


.[. 


° ° ° ( 3 3 



(0) 

(!) (J) 

( 2 4 ) ( 2 5 ) 

(J) (3) 



(% 2 )] 

("T 2 )] 
("")] 
("I 2 )] 



We have E < eP so that E k < e k P k . Let us take the sum of the various powers position by position; 
i.e., find the sum for the (1, 1) position, (1, 2) position, . . ., (1, n) position. 



1,1): 1 

■•MS, 1 



i)- + G) 



{ l^ e+ ' 



1 0) •+G)"+G)" + (a) 



= 1 

= e(l + e)° 

= e(l-he) 2 
= e(l-he) 3 



'-»("o , )- + ("7 , )' + ("7 , )- + --- + (::J)--« + ->' 
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The norm, therefore, comes from adding the extreme righthand column: 

N„[(I -£)-'] = l+e(l + e)° + e(l + c)' + e(l + e)H . . . + e(l + e) w " 2 
= 1 -he[l + (1 -he) + (l + e) 2 + . . . + (1 + e) n ~ 2 ] 

1 -h e 



= (1 -he)"" 1 . 



H 



ence we nave 



= (1 -f e)*" 1 
min \tii\ 

i 

and all the quantities are readily available by inspection of T. 
Corollary 1. 



N„(T- 



min t«t -h max \h*\ 



min |t,i| n 



PROOF: We use the fact that 



tu 



max \ti 



min /,•/ 



for ally 



and substitute this inequality in the bound of the theorem. 

We see in this last bound that the two important numbers are the smallest and the largest. This 
ties in very well with the concept of the condition of a matrix, because a matrix with a large spread 
in the magnitude of its elements is generally considered an ill-conditioned matrix. The larger the 
difference of the maximum element and the minimum element, the higher the value of this bound. 

Corollary 2. 

2 n-i 



//e*l f PUT" 1 ) 



min tii 



This corollary would definitely be applicable to a diagonally dominant triangular matrix. 

Theorem 5.5. Let A be any matrix and let e be defined as in Theorem 5.4. Then if 
N rs (I-AX) «q, 



N/(A-»-X) =S nq 



(1 + e)"- 



min 1 1; | 
Proof: If U is a unitary matrix and T triangular, we have 

A = U*TU 

X = U*RU, for some R 
A- l -X=U*(T- l -R)U 
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sothat N f (A-i-X) = Nf(T-i-R) 

= N f {T~'(I-TR)) 
^NjiT-^Nfil-TR) 

= N s {T-i)N f {I-AX) 

following the same arguments that were used in the proof of Theorem 5.2. 

Now for any A we have the relation N/(A) ^ n li2 N rs {A). Using this inequality we continue 

N f (A-i-X) ^nN rs (T-i)N rs (I-AX). 
Using our bounds for NrsiT' 1 ) we have 

N f (A-i-X)^n K . e \ n Nrsd-AX) 
mm \li\ 

where e= max — for all ;, and U are the eigenvalues of A. The difficulty here, of course, is our 

i tii 

lack of knowledge of ty but the bound certainly brings out the fact that the smaller the eigenvalue 
of A, the higher the bound. This also implies that for a nearly singular matrix (a zero eigenvalue), 
the bounds are much higher which gives another indication of the influence of ill conditioning on 
the calculation of the inverse. 

A more practical application of the bounds for the inverse of a triangular matrix might be made 
to the error in the solution, z, of a linear algebraic system, Ax = b. We pointed out in section 2 that 
to get a bound on the error,^ _1 6 — z, we must use an approximate inverse. This necessitates solving 
an additional n systems. However, if we use the LU decomposition method (Gaussian Elimination), 
we need not do all that. In solving our one given system, the first step is to decompose A into the 
product of an upper triangular matrix t/, and a lower triangular matrix L, whose diagonal elements 
are all ones. 

A = LU 

Then A'^U^L^ 

and N(A-i)*N{U-*)N{L-*). 

Using the bounds just derived for a triangular matrix, we have 

and 



#(£/-*) ^ V T , , , ei=max 

mm \uii\ t 



for ally 



N(L-!)^ K \ 2) , e 2 = max |/ y |. 

1 ij 

The quantities in the last wo bounds are readily available from the decomposition of A and thus 
with very little extra work we will have a bound for the error in the solution vector z without needing 
an approximate inverse. Since for any norm, N(A~ 1 b — z) ^ A^(^ _1 )A^(6~ Az) , we have 
Theorem 5.6. 

N rs (A-ib-z) ^ [(l + e,)(l + e,)r> N r8 (b-Az). 

mm |u«| 

How much of an overestimate this may be again results primarily from using 

N(AB)^N(A)N(B) 
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as mentioned in section 1. We know that the residual is usually much smaller than the actual 
error, however, and this fact may help make this bound practical. 

Examining in itself the LU decomposition of A will give us a slight refinement in the bound 
for L _1 . In the decomposition process we multiply A on the left by a succession of elementary row 
operators to reduce it to triangular form: 



Hence, U- l M n -\M n -2 
so that 

and 



Mn-lMn-2 . . . M t A = U. 
MxA = U- l U = I and therefore 
A-^U-Wn^Mn.2 . . . Mt 

L-^ = M n . l M n -2 ... Mi 
= M 



Now each Mh is a lower triangular matrix with elements in only the A:th column and ones along 
the diagonal: 



M k = 



1 

— mk+\,k 1 



TUnk 1 

so Mk 1 is actually Mu with the signs of the off-diagonal elements reversed. Furthermore 

1 

m 2 i 1 

m u ra 3 2 

m^i TTI42 
L = 



mm m n2 



mn,n-\ 1 
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Now N(Mk) ^ 1 "+■ max \muc\ ^ 1 + max |my| for each A: so that 

i ii 

AT(L-i) =N(M) 



M +max |ra y -| V- 1 



and this is exactly what we have found from our other approach for a triangular matrix. (Note that 
hj — rriij and these are stored in the memory of the computer.) This could be a large overestimate 
since we are repeatedly using the inequality we have indicated as a source of trouble: 



N(M) ^N(M n -i)N(M n - 2 Mn-3 ... MO 



N(M n -l)N(M n -2) • • • N(M t ). 



The refinement comes from not using the second inequality indicated for N(Mk) . Letpj = max 
| ray | for eachy. Then 



and 



Hence we have 



N(M k ) = l+p k 
N(L-i)=N(M) * flf (1+ Pj ). 

N(A~ l ) ^N(L-i)N(U-i) 

* IT (i + pj) r — . 

11 v yj/ nun \uu\ ' 



where e2 — max 

i 

Finally, 



u h 



for ally. 



N(A-*b - z) ^ Yl (1 + B ) A -N(b- 



min \uu\ 



where the notation is the same as before and N is the maximum row sum norm. We therefore have 
a bound on the error of the approximate solution of a linear algebraic system without the knowledge 
of the inverse of the system. 

Let us now give a numerical example of this bound and also of those from section 2. We will 
use the example giver in J.>H. Wilkinson's book Rounding Errors in Algebraic Processes, page 111, 
and we refer the reader to?the discussion given there concerning the size of the bounds for the two 
different values of the righthand sides. 

Problem 1: 



N(A) = .2S2E + 0l 
N(A- 1 ) = A50E + 06 



N(bt) = .912£-h00 N(A-*bi) = .225£ + 05 

N(bi -Az) = .711£-01 N{A~*bx -z) = .151£ + 04 



1. N(A-*)*f[(l+ Pj ) 



{\ + e 2 ) n - 1 
min \ua\ 



[(l + e ,)(l + e 2 )]»- 
min \ua\ 

290 



= P=.861E + 06 

= .112£ + 07 

N(X) = .161E + 06 



2. N(A-*bi-z) *NiA-*)N(bi-Az) =D = .W7E + 05 

orPN(bi-Az) =E=M6E+05 

N(X)N(b,-Az) _ , w+0 , 

»r f_/V(F) -F-.138E + 05 

/V(^)/V(F)/V(^,) _r- W7F-MK 

*(/«-/„ -z) D = 

N(A^bi) N(A-ib t ) -^fi + uu 



£ 



NiA-ifo) 



.274E + 01 



or .,, f ., , = .614E + 00 



iv(i<-*6i; 

G 



A37E + 01 



N(A->b,) 
oiN(I-XA) = .125£+00 
= actually =.673£-01 

4. N(X-*bi-x) ^ N( trJ Z) = .283£-01 
Problem 2: ( ' 

N(A) = .252E+0l N(b-,) = .876E + 00 N(A- l b>) = ,915£+00 

/V(//- 1 ) = .150£ + 06 N(b 2 -Az) = .231E-06 N(A- 1 b-,-z) = .607E-Q2 

(Same as previous problem except that the third element of b\ was decreased by .264139.) 

1. N(A-*b»-z) ^N^A-^Nibi-Az) = £> = .348£-0I 

or PN(bi-Az) =E = .200E+Oi) 

N(X)N(b t -Az) 



l-N(Y) 



= F = .449£-01 



or _\r(Y\ = G = .295L + 05 

^mbt-Az) = 919£ _ 07 

N{A-*b t -z) ^ D _ 

2 - N(A-ibz) *N(A-ib,) -- 380£ - 01 



E 



' N(A-*b t ) 



:.219£ + 00 



M NiA^bT) = - 322 *+ 05 
or/V(/-X4) = .125£+00 
= actually =.668£-02 
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6. Evaluation of Computer Programs 

In this section, the second part of our research, we would like to use the mathematics we 
have developed in our first five sections in the evaluation of the performance of some computer 
programs. We shall set up a problem, choose some programs, and select some test cases to use as 
the data for each program. 

Since calculating the inverse of a matrix usually includes the solution of a linear algebraic 
system, we have chosen the inversion of a matrix for our problem. The programs we will use were 
chosen primarily on the basis of availability. Some were on the computer library tapes at the 
University of Maryland, some were in the literature, and others were received from individuals 
in the field. It seems that almost all of the programs being used today employ some form of Gauss 
elimination to do matrix inversion, so we have confined ourselves to this one direct method. The 
test matrices were chosen to cover a variety of condition numbers, from a small value where all 
programs should perform quite well, to higher values where some difficulties could be encountered. 

Our main concern will be with the accuracy of the result; i.e., with the smallness of N(A~ X — X), 
where X is the approximate inverse produced by a computer program. (This, of course, is not the 
only criterion of success but only the originator of the problem can decide exactly what is important 
to him.) Hence, we use the error bound derived in section 1, that 

7V(/ *" 1-J) " 1 ^{Yy Y = I ~ AX i 

as our basis for evaluation. We will also limit ourselves to the righthand residual matrix since it 
is usually smaller, although, as will be seen, this is not always true. It has been frequently pointed 
out that theoretical error estimates are generally overestimates and sometimes of little practical 
value. For comparison purposes, however, they serve quite well and furthermore, we will show that 
this particular error bound is actually quite close to reality. It is true that we have considered only 
a sampling of computer programs and test problems but we feel they are quite representative and 
that further work will only support our conclusions. 

There are many ways of arriving at a "condition number" for a matrix. The one we shall use 
is called the P- condition number: 



P(A) 



where A is an eigenvalue of largest modulus and /x, of smallest [14] . 

Before going into the details of the test matrices and the computer programs, there are some 
very important items that should be noted first. As we have indicated in many of the previous 
sections, the possibility of the occurrence of some kind of error must always be kept in mind. 
Checking the input, execution, and the output of a computer run may be tedious but it is quite essen- 
tial. We have chosen matrices that can be generated within the machine, used two kinds of checks 
throughout the execution, and have the exact inverses available to check the final results. The main 
check on the execution of the program (and the machine) is achieved by including an extra row and 
column in the original matrix before entering the inversion routine. If the additional column is made 
up of the negative of the sum of the elements of each of the rows of the matrix, then the row sums 
of the augmented matrix should always be zero. If we put a 1 in the (n -h 1, n + 1) position and zeros 
in the other positions of the n 4- 1 row, the last column in the computed inverse should be all ones. 
Let b= (1, 1, . . ., I) 7 . Then 



•Ah]' 1 [A- 1 b 

i J Lo i 



A -AbY 1 [A- 1 b 
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Hence, the maximum difference from unity will give us a good indication if anything has gone wrong. 
What may have happened (machine failure, large rounding errors, etc.) will vary but at least we 
have something definite to tell us to beware of the results. Another check we used was simply 
including an additional column of ones in the second matrix of any product, AB; the result will 
include the row sums of A, which can easily be checked by hand. 

A final check on any run of a computer program should always include a test problem at the 
beginning so that when we examine the output, we can see if the machine performed properly and 
if the program is doing what it is supposed to do. It can also serve as a final check on ourselves; 
for example, if we had dimensioned the variables properly according to all the subroutines we may 
be using. Another check on the output of a computer run comes from a relation that exists between 
the norms we mentioned in the Introduction, the Frobenius norm, Nf, and the maximum element 
norm, Nm> It can be shown that Nf ^ Nm* Hence, using this relation can give us more information 
about the validity of the results. We have found all these checks both necessary and effective. 

One final note on avoiding input errors. In ill-conditioned problems, if the data are not exact, 
the results could be very far from correct. Hence, all our data are integers, exactly represented in 
the machine. This is accomplished by multiplying each matrix by an appropriate scalar, which 
will be indicated in the next subsection. 

The Test Matrices 

This is the fourth power of the 10 X 10 tridiagonal matrix with —2 on the diagonal, 1 above 
and below the diagonal, and elsewhere. 



T=(t li ),t ij = 



h\i~j\ = l 
I 0, \i-j\»2. 



p(T} )~ Wv*y -u 2 ) 4 

*.27£+07. 



The elements are all integers. 

T-i=(b tj ),b tj =i 



-i(n-j+l) 




n + l 


i ^J 


bji 


i >j. 



T is the same as above and we just change the dimension to 20 and the power to 3. 

P(Th) ~ (4/tt 2 ) 3 • (»»)» 
=-.43£ + 07. 
The elements are all integers. 



3. r. 4 



T is the same as above and we just change the power to 4. 



P(Tio) - (4/tT 2 )M" 2 ) 4 

« .69E + 09. 
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The elements are all integers. 

4. A 100 

This is a 10 X 10 matrix where A k = (l/k)I + J and J is the 10 X 10 matrix of all ones. P(A k ) 
= 1 + 10k. The integer form for use as a test matrix is obviously achieved upon multiplication by A:. 

{kA ^ 1=I -iTm J - 

P(A 100 ) = 1001 

« .10£ + 04. 

5. A iooo 

The same as above except that we change the value of k. 

P(A l0 oo) = 10001 

- .10£+05. 

6. A ioooo 

The same as above except that we change the value of k. 

P(A lwm ) = 100001 

- .10£ + 06. 

7. He 

H n is the n X n Hilbert matrix. 

H„=(h ij ),h ij =(i+j-l)-* 
and 

log, [P(H n )] « 3.5/i. 

The integer form for use as a test matrix is achieved upon multiplication by the least common 
multiple of (1, 2, 3, . . ., 2n - 1). 

Hn> = (bij) 

where 

(-l)^(yi + i- l)!(yi+ /-!)! 



ba — 



lJ (i+j-l)[(i-l)l(j-l)l]Hn-i)l(n-j)\ 

P(H 6 ) - .13£+10 

and the multiplier for integer input is 27720. 

8. Hs 
This is the 8x8 Hilbert matrix. 

P(Hs) « .14£+ 13 

and the scalar multiplier for integer input is 360360. 
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9. Hio 

This is the 10 X 10 Hilbert matrix. 

P(#io) «.16£+16 

and the scalar multiplier for integer input is 232792560. 

The Computer Programs 
1. LEQ 

A FORTRAN subroutine used to solve the matrix equation AX = B and to evaluate the determi- 
nant of A. It was written by Max Goldstein of the AEC Computing and Applied Mathematics Center 
at the Courant Institute of Mathematical Sciences, New York University. The Gauss elimination 
method is used. The matrices are normalized row-wise by dividing by the largest element oiA (/, J) 
in that row, then the A matrix is reduced to triangular form by (N — 1) transformations using a 
pivotal condensation process after which X(I, K) is computed by a back-substitution process. This 
transforms B into X and leaves the product of the diagonal elements as the determinant of A. 

2. MATH PACK, GJR 

A FORTRAN subroutine, one of the Univac 1108 Math Pack programs, available on the library 
tapes at the University of Maryland computing center. It solves simultaneous equations, computes 
a determinant, or inverts a matrix or any combination of the three above by using a Gauss-Jordan 
elimination technique with column pivoting. 

3. MATH PACK, MXHOI 

A FORTRAN subroutine, one of the Univac 1108 Math Pack programs, available on the 
library tapes at the University of Maryland computing center. It improves the accuracy of the 
inverse of a matrix by an iterative method which has cubic convergence. 

4. MIDAS 

A FORTRAN and ALGOL package to solve general nonsingular systems of linear algebraic 
equations, invert matrices, and compute determinants. Error bounds on the solution or inverse 
are available as an option. It was written by Peter A. Businger of Bell Telephone Laboratories, 
Inc., Murray Hill, New Jersey. The error bound is a bound on the distance between any element 
of the true inverse and the corresponding element of the computed inverse (unless the bound equals 
— 1, in which case no bound is available). Gaussian elimination with partial pivoting is used to de- 
compose the NxN input matrix into the product of a lower and an upper triangular matrix (LU 
decomposition). The magnitude of intermediate results is estimated; in case of alarming growth 
the program switches to complete pivoting. When solving a system of equations Ax = 6, the accuracy 
of the solution obtained from the triangular systems is improved by iteration; in the case of matrix 
inversion the iteration is omitted for the sake of computational efficiency. (We note that the error 
bound is essentially N(X)N(Y)I[1-N(Y)], Y=I-AX.) 

5. MINV 

A FORTRAN subroutine, one of the IBM System/360 Scientific Subroutine Package. It in- 
verts a general matrix by the standard Gauss-Jordan method. The determinant is also calculated. 
A determinant of zero indicates a singular matrix. 

6. OMNITAB, INVERT 

OMNITAB is a general-purpose computer program for statistical and numerical analysis 
developed at the National Bureau of Standards by Joseph Hilsenrath et al. Now available in an 
ASA FORTRAN version, OMNITAB allows the user to communicate with a computer in an 
efficient manner by means of simple English sentences. INVERT is the calling word for the matrix 
inversion routine used and an error bound is given: N (X)N (Y) I [1 — N(Y)~\. 
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7. SOLVE 

A FORTRAN program by Cleve Moler given in the book Computer Solution of Linear Algebraic 
Systems by George Forsythe and Cleve Moler. It uses Gauss elimination with partial pivoting and 
has a subroutine IMPRUV, which can be called to improve the solution of a linear algebraic system. 
Appropriate messages for various kinds of singularity are available. It is presently undergoing some 
changes to increase efficiency in most FORTRAN systems although these changes should not 
materially alter the numerical behavior. 

8. LSD 

This is a program, with automatic error bounding, for solving a system of linear equations. 
This is equivalent to finding the solution, with error bounds, of the matrix equation Ax= b where 
A is blyi NxN matrix and b and x are N X 1 vectors. The error bounding is done with interval arith- 
metic (I.A.). The program permits interval input and hence errors may be included which arise 
from uncertainty in the data. A second entry to LSD is provided for solving several sets of equa- 
tions with the same A. Thus the program is efficient in finding the inverse of A, with error bounds, 
by using successively for b the columns of the identity matrix. Output from the program includes 
an upper and lower bound on the solution vector. Therefore, the number of significant figures in 
the answer may be simply read from the output. LSD uses two FORTRAN subroutines, a FAP 
function, and a FAP subroutine in addition to interval arithmetic subroutines coded in FAP. The 
program requires only standard FORTRAN II software. The double precision arithmetic hardware 
feature of the IBM 7094 is used, thus restricting the use of LSD, in its present form to this machine. 
It was written by Eldon Hansen and Roberta Smith. 

We intended to include four other programs in our evaluation: 

1. STAT PACK, JIM 

A FORTRAN subroutine, one of the Univac 1108 STAT-PACK programs, available on the 
library tape at the University of Maryland computing center. Apparently, this subroutine never 
functioned properly and the subroutine actually calls GJR of MATH PACK. 

2. SPVMTX 

A single precision FORTRAN IV program for inverting a matrix or solving a set of linear 
equations. To a program from the SHARE library (7090-F1 3180INV1 Single Precision Matrix 
Inversion with Selective Pivoting, written by A. R. Sadaka), Sally T. Peavy, National Bureau of 
Standards, incorporated accuracy checks. We belatedly discovered that this is the routine used 
by INVERT of OMNITAB. This was confirmed by the identical outputs of both programs. 

3. LEQU 

A FORTRAN subroutine to approximate the solution X of the systemAX = B of linear algebraic 
equations, where A is an NX N matrix and B and X are NX M matrices. This program is faster but 
less accurate than subroutine LEQUN, which should normally be preferred. It was written by 
William Kahan, who was then at the Institute of Computer Science, University of Toronto (pres- 
ently at the University of California, Berkeley). We found that this program depends heavily on 
other routines in the University of Toronto 7094 library and we were unable to delete this de- 
pendence or acquire the other routines at this time. 

4. ORTHOl 

An ALGOL subroutine (with iterative improvement) to invert a matrix. It was written by F. L. 
Bauer and essentially uses Gauss elimination but with a suitably weighted combination of rows 
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used for elimination instead of a single row. We were unable, at this time, to run the program in 
ALGOL or to translate it into FORTRAN. 

We would like to mention one other program that exists, one which solves equations exactly. 
It is called SOLVER and was written in FORTRAN by Dr. Morris Newman of the National Bureau 
of Standards. It employs a congruential method and calculates the exact solution of Ax — b where 
the elements of A and b must be integers. It will solve systems in which A is a square matrix at 
most 100 X 100 and the elements of A and b are numerically less than 10 20 . This method is not at 
all sensitive to the condition of A , but it can be time-consuming for large systems. It has successfully 
inverted the Hilbert matrices of all orders up to and including n=\S (a limit imposed by con- 
siderations of time). Since our test matrices all have known inverses, it was not necessary to run 
this program. 

The Results 

In the following tables we list some of the information that we had available from the execution 
of each program for the test problems. We rank the programs according to the smallness of the error 
bound. In some cases, the values of this bound were quite close so the difference in the rank may 
sometimes be slightly artificial. We did not include MXHOI in the rankings and this will be ex- 
plained later. For our estimates we used the Frobenius norm and the maximum element norm and 
we group our results accordingly. In the cases where N (I — AX) was greater than 1, our error bounds 
do not apply and so we leave that row blank. With regard to the check vector carried throughout 
the inversion process, we subtracted this from the vector b= (1, 1, . . ., 1) T and list the maximum 
absolute value of this difference. We note that for a check matrix (T% with a P-condition number of 
approximately 15) all the programs performed extremely well, the only difference, if any, being 
in the last digit. We will discuss the interval arithmetic program of Hansen and Smith separately, 
immediately after presenting the results of the other programs. 

All the calculations, except those for the interval arithmetic program, were performed on the 
Univac 1108 at the University of Maryland between March 10 and March 22, 1970. In use at this 
time was version 23.65.74:N of the EXEC 8 system. 



Table 1. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 

l-N(Y) 


N(X)N(Y) 
l-V(F) 


Relative 
error 


N{I-XA) 


Rank 


LEQ 


.21 £-02 


.13£-01 


.29 £ + 02 


.30 £+03 


.12£-02 


.23 £ + 01 


2 


MATH PAK 


GJR 


.16£-01 


.69 £+00 


.12£ + 04 


.51 £+05 


.86 £-01 


.34 £-01 


6 


MXHOI 


.21 £-02 


.25 £-01 


.66 £ + 02 


.60 £+03 


.29 £-02 


.35 £ + 04 




MIDAS 


.83 £-02 


.29 £-01 


.15£ + 03 


.68 £ + 03 


.65 £-02 


.12 £ + 00 


4 


MINV 


.87 £-02 


.29 £-01 


.15£ + 03 


.69 £ + 03 


.66 £-02 


.27 £-01 


4 


OMNITAB (INVERT) 


.12£-01 


.12£+01 








.46 £-01 


7 


SOLVE 


no IMPRUV 


.35 £-02 


.86 £-02 


.76 £ + 02 


.20 £ + 03 


.33 £-02 


.73 £ + 00 


3 


IMPRUV 


.75 £-08 


.81 £-02 


.15£-03 


.19 £+03 


.65 £-08 


.94 £-02 


1 



Matrix T 4 W 




Condition number 


.27 £ + 07 


Norm Frobenius 



*Difference from 1. 



X= Approximate inverse 



Y=I-AX 



N{A- l -X)-- 
N(A-i-X) 



N(XY) ^N(XJW) 



l-JV(F) l-N(Y) 
NjXY) 1 + /V(F) 



N(A-i) l-N(Y) 



N(X) 
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Table 2. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.20 £-02 


.13£-01 


.22 £ + 03 


.12 £+04 


.25 £-02 


.93 £ + 01 


3 


MATH PAK 


GJR 


.12 £"-02 


.66 £+01 








.20 £-01 


6 


MXHOI 


.41 £-02 


.31 £-01 


.31 £ + 03 


.28 £+04 


.35 £-02 


.65 £ + 04 




MIDAS 


.38 £-02 


.20 £-01 


.24 £ + 03 


.18 £ + 04 


.27 £-02 


.57 £ + 01 


4 


MINV 


.45 £-02 


.31 £-01 


.36 £ + 03 


.29 £ + 04 


.41 £-02 


.23 £-01 


5 


OMNITAB (INVERT) 


.17 £-02 


.70£ + 01 








.36 £-01 


6 


SOLVE 


no IMPRUV 


.17 £-02 


.15 £—01 


.18 £ + 03 


.14 £ + 04 


.21 £-02 


.58 £ + 01 


2 


IMPRUV 


.75 £-08 


.99 £-02 


.60 £-03 


.90 £ + 03 


.68 £-08 


.10£-01 


1 



Matrix T\ 




Condition number 


.43 £ + 07 


Norm Frobenius 



X— Approximate inverse 



Y=I-AX 



N(A~ 
N(A- 



-X) 



N(XY) N(X)N{Y) 



*Difference from 1. 



l-N(Y) l-N(Y) 

X) < N(XY) 1+/V(F) 
N(X) 



N{A^) l-N(Y) 



Table 3. Summary of results 



Program 


Max ELT* 

of check 
vector 


N(I-AX) 


N(XY) 
l-V(F) 


N(X)N(Y) 

\-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.47 £ + 01 


.78 £ + 01 








.22 £ + 05 




MATH PAK 


GJR 


.71 £+00 


.36 £+04 








.31 £ + 01 




MXHOI 




overflow 












MIDAS 


.51 £ + 00 


.28 £ + 01 








.15 £ + 04 




MINV 


.69 £+00 


.32£ + 01 








.21 £ + 01 




OMNITAB (INVERT) 


.37 £+00 


.98 £ + 04 








.45 £ + 01 




SOLVE 


no IMPRUV 


.71 £+00 


.10£ + 01 








.30 £ + 04 




IMPRUV 


.13£-02 


.18 £ + 01 








.14 £ + 02 





Matrix T\ Q 




Condition number 


.69 £ + 09 


Norm Frobenius 



X— Approximate inverse 



Y=I-AX 



N(A-*-X) 



N(XY) N(X)N(Y) 



l-N(Y) l-N(Y) 

N(A-i-X) ^ N(XY) \±NiX) 



*Difference from ] 



N(A~ 



\-N{Y) 



N(X) 
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Table 4. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(f-AX) 


N(XY) 
\-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.75 £-06 


.17 £-04 


.12 £-05 


.50 £-04 


.39 £-06 


.20 £-04 


7 


MATH PAK 


GJR 


.26 £-05 


.42 £-04 


.84 £-06 


.13 £-03 


.28 £-06 


.41 £-04 


6 


MXHOI 


.26 £-05 


.42 £-04 


.84 £-06 


.13 £-03 


.28 £-06 


.41 £-04 




MIDAS 


.86 £-05 


.66 £-04 


.13 £-06 


.20 £-03 


.43 £-07 


.20 £-04 


3 


MINV 


.30 £-07 


.36 £-04 


.39 £-07 


.11 £-03 


.13 £-07 


.19 £-04 


2 


OMNITAB (INVERT) 


.25 £-05 


.41 £-04 


.83 £-06 


.12 £-03 


.28 £-06 


.42 £-04 


5 


SOLVE 


no IMPRUV 


.26 £-05 


.15 £-04 


.82 £-06 


.44 £-04 


.27 £-06 


.18 £-04 


4 


IMPRUV 


.75 £-08 


.72 £-05 


.14 £-07 


.22 £-04 


.47 £-08 


.72 £-05 


1 



Matrix A\w 

Condition number 

Norm Frobenius 



.10 £ + 04 



X= Approximate inverse 



Y=I-AX 



'Difference from 1. 



v ; l-N(Y) l-N(Y) 

N(A-i-X) ^ N{XY) \+N(Y) 
N(A-i) ^l-N(Y) ' N(X) 



Table 5. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N(X)N(Y) 
\-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.34 £-04 


.18 £-03 


.13 £-04 


.54 £-03 


.44 £-05 


.18 £-03 


7 


MATH PAK 


GJR 


.12 £-03 


.33 £-03 


.65 £-05 


.10 £-02 


.22 £-05 


.34 £-03 


4 


MXHOI 


.12 £-03 


.33 £-03 


.65 £-05 


.10 £-02 


.22 £-05 


.34 £-03 




MIDAS 


.67 £-04 


.30 £-03 


.89 £-07 


.91 £-03 


.30 £-07 


.24 £-03 


3 


MINV 


.30 £-07 


.29 £-03 


.36 £-07 


.86 £-03 


.12 £-07 


.18 £-03 


2 


OMNITAB (INVERT) 


.14 £-04 


.43 £-03 


.65 £-05 


.13 £-02 


.22 £-05 


.43 £-03 


4 


SOLVE 


no IMPRUV 


.14 £-04 


.15 £-03 


.65 £-05 


.46 £-03 


.22 £-05 


.16 £-03 


4 


IMPRUV 


.15 £-07 


.84 £-04 


.94 £-08 


.25 £-03 


.31 £-08 


.84 £-04 


1 



Matrix A m w 




Condition number 


.10 £ + 05 


Norm Frobenius 



X= Approximate inverse 



Y=I-AX 



N(A-i-X) 



N(XY) N(X)N(Y) 



'l-N(Y) l-N(Y) 
N(A-i-X) ^ N(XY) \+N(Y) 



*Difference from 1. 



N(A-i) l-N(Y) 



S(X) 



299 









Table 6. Summary of results 








Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 

l-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.42 £"-03 


.13 £-02 


.30 £-03 


.39 £-02 


.10 £-03 


.18 £-02 


7 


MATH PAK 


GJR 


.22 £-06 


.29 J? -02 


.90 £-04 


.88 £-02 


.30 £-04 


.30 £-02 


4 


MXHOI 


.80 £-06 


.54 £-02 


.15 £-03 


.16 £-01 


.51 £-04 


.31 £ + 01 




MIDAS 


.11 £-02 


.27 £-02 


.29 £-06 


.81 £-02 


.96 £-07 


.17 £-02 


3 


MINV 


.30 £-07 


.32 £-02 


.36 £-07 


.95 £-02 


.12 £-07 


.16 £-02 


2 


OMNITAB (INVERT) 


.60 £-07 


.41 £-02 


.90 £-04 


.12 £-01 


.30 £-04 


.41 £-02 


4 


SOLVE 


no IMPRUV 


.30 £-07 


.16 £-02 


.90 £-04 


.48 £-02 


.30 £-04 


.20 £-02 


4 


IMPRUV 


.15 £-07 


.48 £-03 


.10 £-07 


.14 £-02 


.35 £-08 


.48 £-03 


1 



1 10000 



Matrix 

Condition number 

Norm Frobenius 



.10 £ + 06 



A*= Approximate inverse 



Y=I-AX 



N(A-i-X) < 



N(XY) NQCJNJX) 
\-N(Y)* l-N(Y) 

l+N(Y) 



^Difference from 1. 



N(A-'-X) < N(XY) 
N(A-i) 1-N(Y) N(X) 



Table 7. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 

l-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.20 £-01 


.19 £-01 


.92 £ + 00 


.63 £ + 01 


.28 £-02 


.47 £-01 


2 


MATH PAK 


GJR 


.30 £-02 


.30 £-01 


.22 £ + 01 


.10 £ + 02 


.67 £-02 


.31 £-01 


7 


MXHOI 


.11 £-01 


.34 £-01 


.12 £ + 01 


.12 £ + 02 


.36 £-02 


.56 £ + 03 




MIDAS 


.10 £-01 


.28 £ + 00 


.18 £ + 01 


.13 £ + 03 


.67 £-02 


.41 £-01 


6 


MINV 


.86 £-02 


.56 £-01 


.11 £ + 01 


.20 £ + 02 


.36 £-02 


.56 £-01 


3 


OMNITAB (INVERT) 


.10 £-01 


.40 £-01 


.17 £ + 01 


.14 £ + 02 


.52 £-02 


.36 £-01 


5 


SOLVE 


no IMPRUV 


.16 £-01 


.18 £-01 


.16 £ + 01 


.61 £ + 01 


.50 £-02 


.65 £-01 


4 


IMPRUV 


.00 


.96 £-02 


.29 £-05 


.32 £ + 01 


.88 £-08 


.96 £-02 


1 



Matrix # fl 




Condition number 


.13 £+10 


Norm Frobenius 



X— Approximate inverse 



Y=I-AX 



N(A-i-X) *£ 



N(XY) N(X)N(Y) 



'l-N(Y) l-N(Y) 
N(A-i-X) N(XY) 1+N(Y) 



*Difference from 1. 



N(A~ 



\-N{Y) 



N(X) 



300 



Table 8. Summary of results 



Program 


Max ELT* 

of check 

vector 


N{l-AX) 


N(XY) 
l-N(Y) 


N(X)N(Y) 
\-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 




overflow 












MATH PAK 


GJR 


.93 £" + 00 


.20 £ + 02 








.89 £ + 01 




MXHOI 




overflow 












MIDAS 


.77 £ + 01 


.16 £ + 03 








.10 £ + 03 




MINV 


.22 £ + 02 


.26 £ + 02 








.23 £ + 02 




OMNITAB (INVERT) 


.67 £ + 00 


.66 £ + 01 








.52 £ + 01 




SOLVE 


no IMPRUV 


.19 £ + 01 


.20 £ + 01 








.26 £ + 02 




IMPRUV 


.16 £ + 02 


.10 £ + 03 








.35 £ + 03 




Matrix H» 
Condition number 
Norm Frober 

*Difference from 1. 


.14 £ + 


13 


X = Approximate inverse N(A~ l — X) - 
Y I IX N(A-i-X). 


N(XY) ^N(X)N(Y) 

^i-yv(y)" \-N{Y) 


ius 




N(XY) 1+N(Y) 
*l-N(Y) N(X) 



Table 9. Summary of results 



Program 


Max ELT* 

of check 
vector 


N(I-AX) 


N(XY) 
\-N(Y) 


N(X)N(Y) 
l-V(F) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 




overflow 












MATH PAK 


GJR 


.79 £ + 01 


.14 £ + 04 








.67 £ + 02 




MXHOI 




overflow 












MIDAS 


.89 £ + 01 


.35 £ + 02 








.96 £ + 01 




MINV 


.18 £ + 02 


.24 £ + 02 








.22 £ + 04 




OMNITAB (INVERT) 


.17 £ + 02 


.76 £ + 02 








.15 £ + 02 




SOLVE 


no IMPRUV 


.27 £ + 02 


.39 £ + 02 








.16 £ + 04 




IMPRUV 


.47 £ + 03 


.59 £ + 03 








.85 £ + 05 




Matrix / 

Condition nui 
Norm Fr 

*Difference froi 


nber 
oben 

n 1. 


.16 £ + 
ius 


16 


X = A 


sproximate ir 
Y=I-AX 


iverse ^ 


{(A-*-X)s 

'(A-'-X) 
N(A~ l ) 


N(XY) N 
*1-N(Y) * 1 

m N(XY) 1 
e l-N(Y) 


-N(Y) 

+yv(K) 

N(X) 



301 









Table 10. Summary of results 








Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.21 £-02 


.44 £-01 


.53 £ + 02 


.19 £ + 04 


.13 £-02 


.66 £ + 01 


5 


MATH PAK 


GJR 


.16 £"-01 


.18 £ + 01 








.11 £+00 


6 


MXHOI 


.21 £-02 


.91 £-01 


.22 £ + 03 


.42 £ + 04 


.58 £-02 


.86 £+04 




MIDAS 


.83 £-02 


.13 £ + 00 


.29 £ + 03 


.59 £ + 04 


.79 £-02 


.47 £+00 


4 


MINV 


.87 £-02 


.70 £-01 


.27 £ + 03 


.31 £ + 04 


.71 £-02 


.92 £-01 


3 


OMNITAB (INVERT) 


.12 £-01 


.26 £ + 01 








.12 £+00 


6 


SOLVE 


no IMPRUV 


.35 £-02 


.35 £-01 


.14 £ + 03 


.15 £ + 04 


.35 £-02 


.23 £+01 


2 


IMPRUV 


.75 £-08 


.27 £-01 


.46 £-03 


.11 £ + 04 


.11 £-07 


.29 £-01 


1 



Matrix 


7M 

1 10 




Condition 


number 


.27 £ + 07 


Norm 


Maximum el 


ement 



X = Approximate inverse 



Y=I-AX 



N(A~ 
N(A~ 



-X) 



N(XY) N(X)N(Y) 



l-N(Y) i-yv(F) 

X)^ N(XY) 1+V(F) 



*Difference from 1. 



N(A-i) 



l-N(Y) 



N(X) 



Table 11. Summary of results 



Program 


Max ELT* 

of check 
vector 


N(I-AX) 


N(XY) 

l-N(Y) 


N(X)N(Y) 
\-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.20 £-02 


.68 £-01 


.45 £ + 03 


.12 £ + 05 


.28 £-02 


.30 £ + 02 


3 


MATH PAK 


GJR 


.12 £-02 


.23 £ + 02 








.76 £-01 


6 


MXHOI 


.41 £-02 


.13 £ + 00 


.96 £ + 03 


.25 £ + 05 


.63 £-02 


.19 £ + 05 




MIDAS 


.38 £-02 


.77 £-01 


.48 £ + 03 


.14 £ + 05 


.30 £-02 


.28 £ + 02 


4 


MINV 


.45 £-02 


.12 £ + 00 


.75 £ + 03 


.24 £ + 05 


.50 £-02 


.11 £+00 


5 


OMNITAB (INVERT) 


.17 £-02 


.17 £ + 02 








.12 £+00 


6 


SOLVE 


no IMPRUV 


.17 £-02 


.68 £-01 


.37 £ + 03 


.12 £ + 05 


.23 £-02 


.19 £+02 


2 


IMPRUV 


.75 £-08 


.42 £-01 


.22 £-02 


.74 £ + 04 


.13 £-07 


.42 £-01 


1 



Matrix 


n, 




Condition 


number 


.43 £ + 07 


Norm 


Maximum 


element 



X— Approximate inverse 



Y=I-AX 



N(A- l -X) 



N(XY) N(X)N(Y) 



l-N(Y) l-N{Y) 



N(A-i-X) N(XY) 1+/V(F) 



* Difference from ] 



N(A-*) l-N{Y) 



N{X) 



302 









Table 12. Summary of results 








Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
\-N(Y) 


N(X)N(Y) 
l-V(F) 


Relative 
error 


N(I-XA) 


Rank 


LEG 


.47 £ + 01 


.31 £ + 02 








.81 £ + 05 




MATH PAK 


GJR 


.71 £" + 00 


.11 £ + 05 








.11 £ + 02 




MXHOI 




overflow 












MIDAS 


.51 £ + 00 


.12 £ + 02 








.73 £ + 04 




MINV 


.69 £ + 00 


.83 £ + 01 








.79 £ + 01 




OMNITAB (INVERT) 


.37 £ + 00 


.30 £ + 05 








.16 £ + 02 




SOLVE 


no IMPRUV 


.71 £ + 00 


.38 £ + 01 








.15 £ + 05 




IMPRUV 


.13 £-02 


.73 £ + 01 








.44 £ + 02 





Matrix 


1 20 




Condition 


number 


.69 £ + 09 


Norm 


Maximum 


element 



X— Approximate inverse 
Y=I-AX 



N(A-*-X) 



N(XY) N{X)N(Y) 



\-N(Y) \-N(Y) 

N(A- l -X) N(XY) 1 + /V(F) 



*Difference from 1. 



N(A-i) \-N(Y) 



N(X) 









Table 13. Summary of 


results 








Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
\-N(Y) 


N(X)N(Y) 
\-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.73 £-06 


.29 £-04 


.10 £-04 


.26 £-03 


.11 £-05 


.38 £-04 


7 


MATH PAK 


GJR 


.26 £-05 


.69 £-04 


.76 £-05 


.62 £-03 


.84 £-06 


.66 £-04 


6 


MXHOI 


.26 £-05 


.69 £-04 


.76 £-05 


.62 £-03 


.84 £-06 


.66 £-04 




MIDAS 


.86 £-05 


.17 £-03 


.41 £-06 


.16 £-02 


.46 £-07 


.33 £-04 


3 


MINV 


.30 £-07 


.11 £-03 


.12 £-06 


.98 £-03 


.13 £-07 


.20 £-04 


2 


OMNITAB (INVERT) 


.25 £-05 


.57 £-04 


.74 £-05 


.51 £-03 


.82 £-06 


.58 £-04 


5 


SOLVE 


no IMPRUV 


.26 £-05 


.31 £-04 


.73 £-05 


.28 £-03 


.82 £-06 


.25 £-04 


4 


IMPRUV 


.75 £-08 


.76 £-05 


.42 £-07 


.69 £-04 


.47 £-08 


.76 £-05 


1 



Matrix 


Ann) 




Condition 


number 


.10 £ + 04 


Norm 


Maximum 


element 



X— Approximate inverse 



Y=I-AX 



N(A- 



N(A~ 



N(XY) N(X)N(Y) 



' ^l-N(Y) 1-/V(F) 
-X) N(XY) l + N(Y) 



*Difference from 1. 



N(A-*) 1-/V(F) 



N(X) 



303 



Table 14. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N{X)N(Y) 
1-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.34 £-04 


.34 £-03 


.58 £-04 


.31 £-02 


.64 £-05 


.35 £-03 


7 


MATH PAK 


GJR 


.12 £-03 


.57 £-03 


.59 £-04 


.52 £-02 


.65 £-05 


.68 £-03 


4 


MXHOI 


.12 £-03 


.57 £-03 


.59 £-04 


.52 £-02 


.65 £-05 


.68 £-03 




MIDAS 


.67 £-04 


.80 £-03 


.54 £-06 


.72 £-02 


.60 £-07 


.37 £-03 


3 


MINV 


.30 £-07 


.85 £-03 


.93 £-07 


.77 £-02 


.10 £-07 


.18 £-03 


2 


OMNITAB (INVERT) 


.14 £-04 


.63 £-03 


.58 £-04 


.57 £-02 


.65 £-05 


.59 £-03 


4 


SOLVE 


no IMPRUV 


.14 £-04 


.29 £-03 


.58 £-04 


.26 £-02 


.65 £-05 


.27 £-03 


4 


IMPRUV 


.15 £-07 


.88 £-04 


.18 £-07 


.79 £-03 


.20 £-08 


.88 £-04 


1 



Matrix 


^1000 




Condition number 


.10 £ + 05 


Norm 


Maximum element 



X= Approximate inverse 



Y=I-AX 



y ) l-N(Y) l-N(Y) 

N(A-i-X) ^ N(XY) l + N(Y) 



*Difference from 1. 



N(A~ l ) l-N(Y) 



N(X) 



Table 15. Summary of results 



Program 


Max ELT* 

of check 
vector 


N(I-AX) 


N{XY) 
l-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.42 £-03 


.34 £-02 


.21 £-02 


.31 £-01 


.24 £-03 


.30 £-02 


7 


MATH PAK 


GJR 


.22 £-06 


.59 £-02 


.82 £-03 


.53 £-01 


.91 £-04 


.65 £-02 


5 


MXHOI 


.80 £-06 


.96 £-02 


.70 £-03 


.87 £-01 


.79 £-04 


.63 £ + 01 




MIDAS 


.11 £-02 


.60 £-02 


.12 £-05 


.54 £-01 


.13 £-06 


.25 £-02 


3 


MINV 


.30 £-07 


.98 £-02 


.11 £-06 


.89 £-01 


.12 £-07 


.16 £-02 


2 


OMNITAB (INVERT) 


.60 £-07 


.61 £-02 


.82 £-03 


.55 £-01 


.91 £-04 


.62 £-02 


5 


SOLVE 


no IMPRUV 


.30 £-07 


.29 £-02 


.81 £-03 


.26 £-01 


.91 £-04 


.40 £-02 


4 


IMPRUV 


.15 £-07 


.49 £-03 


.31 £-07 


.45 £-02 


.35 £-08 


.49 £-03 


1 



Matrix 


^10000 




Condition 


number 


.10 £ + 06 


Norm 


Maximum element 



X = Approximate inverse 



Y=I-AX 



N(A-i-X) *s 



*Difference from 1. 



N(XY) _ N(X)N(Y) 
l-N(Y)^ l-N(Y) 

N(A-i-X) < N(XY) 1+N(Y) 
N(A-i) ^l-JV(F) ' N(X) 



304 









Table 16. Summary of results 








Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N(X)N(Y) 
l-N(Y) 


Relative 
error 


N(I-XA) 


Rank 


LEQ 


.20 £"-01 


.57 £-01 


.26 £ + 01 


.58 £ + 02 


.29 £-02 


.12 £ + 00 


2 


MATH PAK 


GJR 


.30 £-02 


.12 £ + 00 


.68 £ + 01 


.14 £ + 03 


.80 £-02 


.97 £-01 


6 


MXHOI 


.11 £-01 


.11 £ + 00 


.32 £ + 01 


.12 £ + 03 


.37 £-02 


.12 £ + 04 




MIDAS 


.10 £-01 


.86 £ + 00 


.27 £ + 02 


.61 £ + 04 


.52 £-01 


.16 £ + 00 


7 


MINV 


.86 £-02 


.19 £ + 00 


.37 £ + 01 


.23 £ + 03 


M £-02 


.22 £ + 00 


3 


OMNITAB (INVERT) 


.10 £-01 


.91 £-01 


.51 £ + 01 


.96 £ + 02 


.58 £-02 


.12 £ + 00 


5 


SOLVE 


no IMPRUV 


.16 £-01 


.61 £-01 


.49 £ + 01 


.62 £ + 02 


.54 £-02 


.13 £ + 00 


4 


IMPRUV 


.00 


.23 £-01 


.11 £-04 


.22 £ + 02 


.11 £-07 


.23 £-01 


1 



Matrix 



H, 



13 £+10 



Condition number _ 

Norm Maximum element 



X = Approximate inverse 



Y=I-AX 



N(A~ l -X) 
N(A-*) 



NjXY) ^.N(X)N(Y) 



l-N(Y) \-N(Y) 



N(XY) 
l-N(Y) 



l + N(Y) 
N(X) 



*Difference from 1. 



Table 17. Summary of results 



Program 


Max ELT* 

of check 
vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N(X)N(Y) 

l-N(Y) 


Relative 
error 


N{I-XA) 


Rank 


LEQ 




overflow 












MATH PAK 


GJR 


.93 £ + 00 


.67 £ + 02 








.25 £ + 02 




MXHOI 




overflow 












MIDAS 


.77 £ + 01 


.49 £ + 03 








.67 £ + 03 




MINV 


.22 £ + 02 


.81 £ + 02 








.80 £ + 02 




OMNITAB (INVERT) 


.67 £ + 00 


.29 £ + 02 








.16 £ + 02 




SOLVE 


no IMPRUV 


.19 £ + 01 


.80 £ + 01 








.66 £ + 02 




IMPRUV 


.16 £ + 02 


.26 £ + 03 








.12 £ + 04 





Matrix 


H 8 




Condition 


number 


.14£+13 


Norm 


Maximum 


element 



X= Approximate inverse 



Y=I-AX 



N(A-i-X) 



N{XY) N(X)N(Y) 



l-N(Y). l-N(Y) 
N(A-i-X) N(XY) 1+/V(F) 



*Difference from 1.' 



A^- 1 ) l-N(Y) 



N(X) 



305 



Table 18. Summary of results 



Program 


Max ELT* 

of check 

vector 


N(I-AX) 


N(XY) 
l-N(Y) 


N{X)N{Y) 
l-N(Y) 


Relative 
error 


• N(I-XA) 


Rank 


LEQ 




overflow 












MATH PAK 


GJR 


.79 £" + 01 


.58 £ + 04 








.26 £ + 03 




MXHOI 




overflow 












MIDAS 


.89 £" + 01 


.12 £ + 03 








.61 £ + 02 




MINV 


.18 £" + 02 


.12 £ + 03 








.94 £" + 02 




OMNITAB (INVERT) 


.17 £ + 02 


.20 £ + 03 








.79 £ + 02 




SOLVE 


no IMPRUV 


.27 £ + 02 


.13 £ + 03 








.40 £ + 04 




IMPRUV 


.47 £ + 03 


.24 £ + 04 








.25 £ + 06 





Matrix 


Hu) 




Condition 


number 


.16 £+16 


Norm 


Maximum 


element 



X = Approximate inverse 



Y=I-AX 



N(A-*-X) : 



N(A~ l -X) 



N(XY) N(X)N(Y) 



l-N(Y) l-N(Y) 



* Difference from 1. 



N(XY) 1+N(Y) 
N(A~ l ) ^l-N(Y) ' N(X) 



For the program LSD, using interval arithmetic, we found excellent results for all the matrices 
except H$, //io, and Too. By excellent results we mean the upper and lower values given for each ele- 
ment of the inverse agreed to 7 or 8 digits with each other and with the exact value. For //g, the 
spread between the upper and lower values was quite different. A typical situation is given in the 
first row below, and the worst spread in the second row. 



A^U 



.12533643 £"-03 
.16101655 £ + 01 



Exact 
.17760000 £-03 
.80000000 E+ 00 



A~ x Upper 
.22986282 £-03 
+ .10338515 £-01 



In view of the fact that none of the other programs solved this problem, these results are not bad. 

For the 10 X 10 Hilbert matrix, the following error message was printed out: SINGULAR 
MATRIX INCLUDED IN INPUT MATRIX. (We note that SOLVER produced the exact inverse 
in both of these cases.) 

For Tic we got results but the spread was so large that they are meaningless. One of the better 
values was from .26055424 £ + 05 to .56997728 £ + 05 but a more common spread would be similar 
to the following: 



or 



-.14807597 £ + 07 to 
-.60322834 £ + 08 to 



+ .18754179 £ + 07 
+ .60828502 £ + 08. 



These certainly are not very satisfactory. Let us now turn our attention to the other programs. 

Evaluation and Discussion 

It is obvious that SOLVE with iterative improvement performed far better than any of the others. 
The excellent results that are achieved by using this added feature certainly seem to warrant the extra 
work involved. In each case where it did not succeed in inverting the matrix (H H , H U) , r| ), the follow- 
ing error message was printed out: ITERATION DID NOT CONVERGE. MATRIX IS NEARLY 
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SINGULAR. This is most important because it tells us our answers may not be correct. LEQ had 
floating point overflow, MIDAS had the error bound of — 1, and OMNITAB had a negative bound, 
each of which also indicates doubtful results. MINV and GJR gave no indication that the results 
might not be correct— a very serious defect. GJR does, in fact, have an error return but it is only for 
overflow and this is not sufficient. The zero determinant criterion of MINV is also not adequate. 
Knowing that the condition number is as high as it is for these two matrices, we expected trouble. 
But for T20 whose condition number is fairly large but whose largest element is 70, we did not antici- 
pate serious difficulty and the message of singularity was needed. In this case, however, LEQ did 
not have any floating point overflow message. 

If we did not know ahead of time that the matrix was poorly conditioned, we could use the 
iterative improvement scheme if for no other reason than to inform us of doubtful results. In solving 
a single system of equations where we may not calculate an approximate inverse (from which a 
residual matrix whose norm being greater than 1 can indicate trouble in spite of the program out- 
put), the improvement scheme seems an excellent way of indicating possible poor results. 

In the cases where the inversions were successfully accomplished (that is, at least N(I — AX) 
was less than 1), it seems that MINV, SOLVE (no IMPRUV), and MIDAS were more successful 
than the others, excluding IMPRUV of course. It is interesting to note that the only program (besides 
IMPRUV) which had the norm of both residuals less than one was MINV. This is certainly a desir- 
able situation. It also points out the large discrepancy that can exist between the two residual 
matrices and the fact that the righthand residual,/ — AX , is not always smaller. Hence, a bound on 
the error matrix can sometimes be available by using the other residual. None of the programs 
seem to take this possibility into consideration. If we have gone to the trouble of calculating/ — AX 
and get a negative bound, it might be worthwhile to calculate/ — XA and see if a bound can be found. 
This will also indicate the kind of inverse we have. It is still, of course, up to the originator of the 
problem to decide if this is necessary. 

The vagaries of numerical calculations are also clearly indicated by these results. We have 
cases in the same problem where N(I — AX) > N(I — XA) for one program and the opposite is true 
for another program; for the same program, we may have N(I — AX) > N(I — XA) for one problem 
and the opposite for another problem. Hence, one must always be extremely careful in doing 
numerical work. 

The poor and erratic performance of MXHOI led to a brief investigation of this improvement 
procedure. It never really improved the results and in some cases (/4ioooo) made them worse. We 
found that no double precision was used at all and, as we indicated in section 3, for a scheme of 
cubic convergence, multiprecision may be needed, let along double precision. The performance of 
MXHOI certainly bears this out! The results are completely unreliable and this is why we did not 
even consider this routine as far as ranking the various programs was concerned. 

Finally, we would like to discuss briefly the merits of the bound N (XY)/ [I— N(Y)]. We have 
observed in section 1 that using this bound would probably give better results for the relative error 
than the simpler expression /V(F), where Y is either residual matrix. We note that this is in fact 
true in every one of our problems where the inversion was successful. 

Another and much more important fact is the following. We have attempted in our first part 
to derive a practical bound for N(A~ l —X); i.e., one that is realistic and not misleading. It is clear 
that using N(XY) is better than N(X)N(Y) —there is almost always at least a factor of 10 involved 
in the improvement. This is certainly worthwhile. But just how close to the actual error is our bound? 
This question can be partially answered by looking at the lower bound, N(E)/2N(A) . This was not 
included in the tables because the only time it was close to the upper bound was for SOLVE with 
IMPRUV, which certainly indicates both the effectiveness of this scheme and the reality of the 
upper bound. A more complete answer comes from calculating the actual error. We chose He and 
T$ as being of the more difficult type and list the results in the following tables. The comparisons 
are extremely satisfying. 
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TABLE 19. Norms of actual difference 









Maximum Element Norm 


Frobenius Norm 


Program 


N(XY) 
l-N(Y) 


N{A-!-X) 
(actual) 


N(X)N(Y) 
l-N(Y) 


N(XY) 
l-N(Y) 


N(A~*-X) 
(actual) 


N(X)N(Y) 
l-N(Y) 


LEQ 


.26 £ + 01 


.25 £ + 01 


.58 £ + 02 


.92 £ + 00 


.91 £ + 00 


.63 £ + 01 


MATH PAK 


GJR 


.68 £ + 01 


.59 £ + 01 


.14 £ + 03 


.22 £ + 01 


.21 £ + 01 


.10 £ + 02 


MXHOI 














MIDAS 


.27 £ + 02 


.37 £ + 01 


.61 £ + 04 


.18 £ + 01 


.13 £ + 01 


.13 £ + 03 


MINV 


.37 £ + 01 


.30 £ + 01 


.23 £ + 03 


.11 £ + 01 


.11 £ + 01 


.36 £ + 02 


OMNITAB (INVERT) 


.51 £ + 01 


.46 £ + 01 


.96 £ + 02 


.17 £ + 01 


.16£ + 01 


.14 £ + 02 


SOLVr 


no IMPRUV 


.49 £ + 01 


.46 £ + 01 


.62 £ + 02 


.16 £ + 01 


.16£+01 


.61 £ + 01 


IMPRUV 


.11 £-04 


.60 £-06 


.22 £ + 02 


.29 £-05 


.14 £-06 


.32 £+01 



Matrix . 



Condition number . 



JU. 



.13 £+10 



X = Approximate inverse N{A~ l — X) 
Y=I-AX 



\-N{Y) 

N(A-'-X)^ N(XY) 



N(XY) ^ N(X)N(Y) 

l-N(Y) 



N(A~ l ) l-N(Y) 



1 + N(Y) 

N(X) 



Table 20. Norms of actual difference 









Maximum Element 


Norm 


Frobenius Norm 


Program 


N(XY) 
l-N(Y) 


N(A-'-X) 
(actual) 


N(X)N(Y) 
l-N(Y) 


N(XY) 
l-N(Y) 


N(A-'-X) 
(actual) 


N(X)N{Y) 
\-N{Y) 


LEQ 


.45 £ + 03 


.42 £+03 


.12 £ + 05 


.22 £ + 03 


.22 £ + 03 


.12 £ + 04 


MATH PAK 


GJR** 


.40 £ + 03 


.37 £+03 


.14 £ + 05 


.19 £ + 03 


.19 £+03 


.18 £ + 04 


MXHOI 














MIDAS 


.48 £ + 03 


.44 £+03 


.14 £ + 05 


.24 £+03 


.23 £+03 


.18£ + 04 


MINV 


.75 £ + 03 


.66 £ + 03 


.24 £ + 05 


.36 £ + 03 


.35 £ + 03 


.29 £ + 04 


OMNITAB (INVERT)** 


.39 £ + 03 


.36 £ + 03 


.24 £ + 05 


.19 £ + 03 


.18 £ + 03 


.34 £ + 04 


SOLVE 


n6 IMPRUV 














IMPRUV 


.22 £-02 


.12 £-02 


.74 £ + 04 


.60 £-03 


.36 £-03 


.90 £ + 03 



Matrix . 



X = Approximate inverse N(A~ 1 —X)' 



NjXY) N(X)N(Y) 
l-N(Y)^ l-N(Y) 



Condition number . 



,43 £ + 07 



Y=I-AX 
*Y=I-XA 
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N(A-*-X) ^ N(XY) l + N(Y) 
N(A~ l ) " \-N(Y) ' N(X) 



In conclusion, we would like to make the following observations: 

1. Some indication of the possibility of poor results should be included in the output of every 
program. 

2. Some indication of the accuracy of the results should be included in the output of every 
program. 

3. SOLVE with iterative improvement is certainly an excellent program. 

4. Interval arithmetic was quite effective but not in every case. 

N(XY) 

5. The use of _^/yx gives a realistic (practical) upper bound on the error of a computed 

inverse, can be used with either residual matrix and would indicate the possibility of poor results 
if the norm of each residual matrix were greater than one. 



I wish to thank the many people who encouraged me throughout my pursuit of this final step, 
especially Morris Newman, who made it all possible. 

I am grateful to the Lewis and Rosa Strauss Memorial Fund for the financial aid which made 
the computing facilities of the University of Maryland available to me. 
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